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1. Introduction 

The POWHEG method is a prescription for interfacing NLO calculations with parton shower 
generators. It was first suggested in ref. [1], and was described in great detail in ref. [2]. 
Until now, the POWHEG method has been applied to the following processes: Z pair hadropro- 
duction [3], heavy-flavour production [4], e + e~ annihilation into hadrons [5] and into top 
pairs [6], Drell-Yan vector boson production [7, 8], W' production [9], Higgs boson pro- 
duction via gluon fusion [10, 11], Higgs boson production associated with a vector boson 
(Higgs-strahlung) [11], single-top production [12], Z + 1 jet production [13], and, very re- 
cently, Higgs production in vector boson fusion [14]. Unlike MC0NL0 [15], POWHEG produces 
events with positive (constant) weight, and, furthermore, does not depend on the Monte 
Carlo program used for subsequent showering. It can be easily interfaced to any modern 
shower generator and, in fact, it has been interfaced to HERWIG [16, 17] and PYTHIA [18] in 
refs. [3, 4, 7, 10, 12, 14]. 

The present work introduces a computer framework that implements in practice the 
theoretical construction of ref. [2]. We call this framework the POWHEG BOX. The aim of 
the POWHEG BOX is to construct a POWHEG implementation of an NLO process, given the 
following ingredients: 
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i. The list of all flavour structures of the Born processes. 

ii. The list of all flavour structures of the real processes. 

iii. The Born phase space. 

iv. The Born squared amplitudes B, the color correlated ones Bij and spin correlated ones 
B^ u . These are common ingredients of NLO calculations performed with a subtraction 
method. 

v. The real matrix elements squared for all relevant partonic processes. 

vi. The finite part of the virtual corrections computed in dimensional regularization or in 
dimensional reduction. 

vii. The Born colour structures in the limit of a large number of colours. 

With the exception of the virtual corrections, all these ingredients are nowadays easily 
obtained. A matrix element program (like MADGRAPH) can be used to obtain (iv) and (v). 
The colour-correlated and spin-correlated Born amplitudes are also generated automatically 
by programs like MadDipole [19] and AutoDipole [20]. Recent progress in the automatiza- 
tion of the virtual cross section calculation may lead to developments where even the virtual 
contribution (vi) may be obtained in a painless way [21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. 
Given the ingredients listed above, the POWHEG BOX does all the rest. It automatically finds 
all the singular regions, builds the soft and collinear counterterms and the soft and collinear 
remnants, and then generates the radiation using the POWHEG Sudakov form factor. 

The purpose of this work is twofold. Our first aim is to complete here the theoretical 
work of ref. [2], by explaining several variants of the procedure illustrated there, that have 
turned out to be necessary to fulfill our goal. In doing so, we will refer often to formulae 
given in ref. [2], that is thus a prerequisite for reading the present work. Our second 
aim is to illustrate specifically how the code really works. It is true to some extent that 
well-written codes are self explanatory, and, in fact, we tried to write the POWHEG BOX as 
transparently as possible. However, what may not be so easy to understand from the code 
is its global organization. We believe that this document, together with the source code, 
could be used by others to understand the code up to the point of being able to modify it. 

Strictly speaking, the present work is neither the documentation of the POWHEG BOX 
code, nor a description of its theoretical basis. So, for example, we do not include a 
rigorous description of all the subroutines in the code, and, as we already said, we refer to 
ref. [2] for the theoretical bases of our method. In general, one expects that a theoretical 
paper should include the description of how a given calculation has been performed. This 
normally includes the illustration of some algebraic steps, and, maybe, an indication of 
what part of the calculation was performed numerically, so that the reader should be 
able, on the basis of the given indications, to verify its content. In the present case, the 
calculation is performed relying heavily on computer algorithms, and the only realistic 
way to verify its content is to understand the code. Thus, here we explain the algorithms 
and give sufficient indications on the code structure for the reader to understand it and 
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verify its correctness. We believe that a detailed documentation of the POWHEG BOX is not 
necessary for this purpose, since the details are better understood by directly studying 
the code. Furthermore, the code itself will unavoidably evolve with time, so that detailed 
documentation may not be so helpful after all. In summary, this paper should be simply 
seen as a description of our calculations, that, being performed essentially by a computer 
program, must, to some extent, coincide with a description of the program itself. 

Researchers wishing to use the POWHEG BOX should not need, in principle, to study or 
understand this whole paper. They only need to know in which format they have to supply 
the ingredients that are listed above. These are summarized in section 2. Looking at the 
various implementation examples included in the code should also help with this task. In 
the near future, we will provide a manual that documents in detail the user interface. The 
remaining part of this work should be useful in order to understand better certain features 
that the POWHEG BOX provides, and that the user may need, or to implement new features 
that are not yet available. 

The paper is organized as follows: in section 2 we report all the information needed 
to interface an NLO program to the POWHEG BOX. Thus, for example, we specify how the 
flavour structure of scattering processes is represented, how to specify their kinematics, 
and the format that the Born, virtual and real cross section subroutines must have. In 
section 3 we describe the algorithm used to generate all the singular regions of the process, 
and how these are represented in the computer code. A simple example is also discussed 
in detail. In section 4 we describe how the B function is constructed. This function, when 
integrated over the radiation variables, yields the inclusive cross section at fixed underlying 
Born configuration, and thus is crucial for the first stage of the event generation, i.e. the 
generation of the underlying Born structure. The computation of B is quite complex, 
so this section is divided into several subsections, dealing with the Born contribution, the 
soft-virtual contribution, the real contribution and its soft and collinear limits. In section 5 
we describe the mechanism provided in POWHEG BOX for tuning the part of the real cross 
section that is dealt with using the shower Monte Carlo method, and how the remaining 
finite part is treated. This separation into a singular and a finite part, besides being useful 
for tuning the POWHEG output, also provides a method to improve generation efficiency in 
certain cases. In section 6 we give an overview of the initialization stage of POWHEG. At 
this stage the inclusive cross section is computed, and the appropriate grids are set up for 
the generation of the underlying Born configurations. In section 7 we describe the second 
stage of the event generation, i.e. the generation of radiation, and in section 8 we describe 
how the POWHEG-generated event is prepared for further showering by a standard shower 
Monte Carlo program. Finally, in section 9, we give our conclusions. 

Several appendixes collect further analytical and technical details. In appendix A we 
report the formulae for the soft integrals that are used in the POWHEG BOX. In appendix B 
we report the collinear limits of the real cross section, that are used in the POWHEG BOX to 
build the collinear counterterms. In appendixes C and D we describe the upper bounding 
functions used in the generation of radiation, both for final-state and for initial-state ra- 
diation. In appendix E we describe how the renormalization and factorization scales are 
set in the POWHEG BOX, and how the strong coupling constant is computed. Finally, in 
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appendix F, we give a discussion of a few miscellaneous topics that are useful for using and 
understanding the program. 

2. The format of the user subroutines 

By flavour structure of a process we mean the type of all incoming and outgoing particles. 
In the program, the flavour type is denoted by an integer number, so that the flavour 
structure is a list of integers. The ordering and numbering of particles follow the rules: 

1. first particle: incoming particle with positive rapidity 

2. second particle: incoming particle with negative rapidity 

3. from the third particle onward: final-state particles ordered as follows 

— colourless particles first, 

— massive coloured particles, 

— massless coloured particles. 

The flavour is taken incoming for the two incoming particles and outgoing for the final-state 
particles. 

The flavour index is assigned according to the Particle Data Group conventions [31], 
except for the gluons, where (rather than 21) is used. 

Example: if we are interested in the associated production of a t quark and a vector boson 
Z plus two jets 

pp^Zt + 2 jets, (2.1) 
then one of its contributing subprocesses is 

bu^Ztsg, (2.2) 

whose flavour structure, according to the previous rules, is given by 

[5, 2, 23, 6, 3, 0] . (2.3) 

In QCD calculations, the colourless particles and the massive coloured particles will remain 
the same at the Born and NLO level. In the NLO calculation, the flavour structure of real 
graphs will have one more light parton in the final state. The virtual term, being the 
interference of the Born and of the one-loop amplitude, has the same flavour structure of 
the Born term. 

In the PDWHEG BOX, the header file pwhg_f 1st .h, in the include subdirectory, contains 
all arrays and parameters having to do with flavour structures. They depend upon the 
parameters nlegborn and nlegreal that are set to the number of legs (incoming plus 
outgoing) of the Born (or virtual) and real graphs. The user of the POWHEG BOX will 
have to set explicitly nlegborn in the nlegborn. h file, that is thus included before the 
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pwhg_f 1st .h file in all program units that need to access the flavour structures. In the 
example (2.2), the user should set nlegborn=6. The variable nlegreal is always set to 
nlegborn+1. The user should also set the variables flst_nborn and flstjireal to the 
number of inequivalent flavour structures for the Born and real graphs, and should also fill 
the arrays 

f lst_born(k=l :nlegborn, j = l :f lst_nborn) 
f lst_real(k=l :nlegreal, j = l : flst_nreal) 

in an appropriate initialization subroutine, named init_processes. Notice that flavour 
structures that are equivalent under a permutation of final-state particles should never 
appear in the list. Thus, in the example of eqs. (2.2) and (2.3), either [5, 2, 23, 6, 3, 0] or 
[5, 2, 23, 6, 0, 3] may appear as flavour structures, but not both. 

The user should also set the value of f lst_lightpart, the position of the first light 
coloured parton. Then, in the example previously described, f lst_lightpart=5. 

2.1 Tagging parton lines 

At times it is convenient to treat lines with the same flavour as if they were different. 

One such example is Higgs boson production via vector boson fusion (VBF). The 
fermion lines attached to the vector bosons may be treated as being distinct. Of course, 
in W + W~ fusion, they may also be effectively distinct, with the W + coming from a u 
quark turning into a d, and the W~ coming from a c quark turning into an s. Consider 
however the real graph depicted in fig. 1. It corresponds to a gluon- initiated next-to-leading 




Figure 1: Example of NLO gluon- initiated correction to Higgs boson production in VBF: gu — > 
Hudd. 

correction to VBF Higgs boson production: gu — > Hudd. It is clear that the two d quarks 
in the final state have a very different role, and should be kept distinct. However, as 
far as the flavour combinatorics is concerned, they are considered identical in the POWHEG 
BOX, that assumes that the graphs are already symmetrized with respect to identical final- 
state particles. Thus, the combinatoric algorithm will generate two regions for this graph, 
corresponding to either d being collinear to the incoming gluon. In order to overcome this 
problem, the POWHEG BOX allows the possibility to attribute a tag to each line, so that lines 
with the same flavour but different tags will be treated differently from the combinatoric 
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point of view. In the example at hand, one assigns the tags according to the scheme in 
fig. 2. We arbitrarily assign a tag equal to zero to particles that we do not need to tag (the 
initial-state gluon and the produced Higgs boson). Within this scheme, the two final-state 
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Figure 2: Tag assignment for the underlying Born graph du — > Hud and its gluon- initiated real 
diagram gu — > Hudd. 

d quarks will be treated as different from the combinatoric point of view. Only the quark 
tagged as 1 in the real graph will generate a singular region, since if quark 2 were collinear 
to the incoming gluon, the associated underlying Born would have an incoming antiquark 
tagged as 2, and thus would not be present. 

In the code, tagged lines are treated by generating an internal flavour index that 
replaces the real flavour index, in such a way that the internal flavour is different for lines 
with different flavour or different tag. The combinatorics is carried out with these internal 
flavour numbers. At the end, internal flavour numbers are replaced with the original real 
flavour numbers. From this point on, tags are ignored. The tags are set to zero during 
POWHEG initialization. If needed, the user's initialization routine should appropriately set 
the arrays f lst_borntags and f lst_realtags. 

The task of changing the flavour indexes to account for different tags, and of chang- 
ing them back to the original state, are carried out by the routines mapflavours and 
unmapf lavours, that are called near the beginning and near the end of the genf lavreglist 
routine, the subroutine that finds the different singular regions, as described in section 3. 

2.2 The Born phase space 

The Born phase space (provided by the user of the PDWHEG BOX) is a subroutine named 
born_phsp (xborn) that fills the four-momenta of Born-process particles (both in the lab- 
oratory and in the center-of-mass frame), the Bjorken x of the two incoming partons, the 
minimal mass of the Born system and the phase space volume. 

It receives as input xborn, an array of real numbers, in the range [0, 1]. If no resonances 
are present in the final state, the dimension of this array is 3n — 2, n being the number 
of final-state particles, i.e. n =nlegborn— 2 . In case some resonances are present, their 
virtuality will require one more variable, and the user will have to take account of them 
too, increasing accordingly the numbers of elements of this array. 
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In the paper, we will refer to this set of numbers as -X"bom- We recall that the Born 
phase space <& n , defined in [2], is given by 



/ n \ n ,3, 

d*„ = dx® dx e (2ir)W fc e + k e - £ h J] l 2k0 . (2.4) 



The born_phsp routine should perform the following tasks: 

1. Set kn_pborn (mu=0 : 3 , k=l :nlegborn) and kn_cmpborn(mu=0 : 3 , k=l :nlegborn) 1 to 
the Born momenta in the laboratory frame and in the center-of-mass (CM) frame. The 
Lorentz index fi = denotes the time component, 1,2 the transverse components x,y, 
and 3 the longitudinal component z. Set the variables kn_xbl and kn_xb2 to the value 
of the parton momentum fraction and x Q . Set the variable kn_sborn to the squared 
CM energy of the Born process. 

2. The array kn_masses should be filled with the masses of the legs of the process. Fur- 
thermore, the variable kn_minmass should be set to a fixed (i.e. independent upon the 
kinematics) lower bound on the mass of the final state. Thus, if no resonances are 
present, it is typically set to the sum of the masses of the final-state particles. If there 
are resonances, it will be set to the sum of the lower limits of the windows imposed 
around the resonances. 



3. Set the variable kn_j acborn to the Jacobian 



Jbom — 



<9^born 

2.3 The Born and Born-correlated squared amplitudes 



(2.5) 



The user of the POWHEG BOX should provide the routine 

setborn(p(0:3, 1 :nlegborn) , bf lav(l :nlegborn) , born, 

bornjk(l :nlegborn, 1 :nlegborn) , bmunu(0 : 3 , : 3 , 1 :nlegborn) ) . 

Given the four-momenta p and the flavour structure bf lav of a Born subprocess, the routine 
should return the Born squared matrix element Is^B in born, the colour correlated one in 
bornjk and the spin correlated one in bmunu. The flux factor 1/(2 Sf,) =1/ (2*kn_sborn) 
(where Sb is the center-of-mass energy squared of the Born process) should not be included, 2 
since it is supplied by the POWHEG BOX. 

The colour correlated Born amplitude is defined in eq. (2.97) of ref. [2]. We report it 
here for completeness 

2s h B i3 =~NJ2 M {Ck} (M\ Ck} ) c ^ c , T c ; c , . (2.6) 

spins c j~* c j 
colours 

Here M{ Ck y is the Born amplitude, and {c^} stands for the colour indexes of all external 
coloured particles in the amplitude. The suffix on the parentheses that enclose -M{ Cfc } 



X A11 variables with the kn_ prefix are defined in the header file pwhgjtn.h. 
2 In the notation of ref [2], B includes the flux factor 
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indicates that the colour indexes of partons i,j are substituted with primed indexes in 
A^| Cfc |. The factor N is the appropriate normalization factor including averages over initial 
spin and colour and symmetry factors. We assume summation over repeated colour indexes 
(cfc, c^, c'j and a) and spin indexes. For gluons T® b = if cab, where f a b c are the structure 
constants of the SU(3) algebra. For incoming quarks = t^g, where t are the colour 
matrices in the fundamental representation (normalized as Tr [tt] = 1/2). For antiquarks 
T£g = —t a 0a . It follows from colour conservation that Bij satisfy 

Y,B ij = C fj B, (2.7) 

where % runs over all coloured particles entering or exiting the process, and Cf i is the 
Casimir constant for the colour representation of particle j. The spin correlated Born 
squared amplitude B^ u is defined to be non-zero if the j th Born leg is a gluon, and is 
basically the Born cross section obtained by leaving the gluon indexes of the j th leg un- 
contr acted. More precisely, we can write 

Bf = N £ M({i},s j )M^({i},s' j ) (e^reZj, (2.8) 

where M. ({i}, Sj) is the Born amplitude, {i} represent collectively all remaining spins and 
colours of the incoming and outgoing particles, and Sj represents the spin of the j th particle. 
The are polarization vectors, normalized as 

Y,%ru>: )■<:■ s SjS ,. (2.9) 

Thus 

Y,9,uB^ = -B. (2.10) 

Notice that the Born squared amplitude is requested for each individual flavour structure 
of the contributing subprocesses. Many different flavour structures will return identical or 
proportional values of the Born cross section. For example dd — > Z is identical to ss — > Z, 
and uu — > 7* is proportional to dd — >■ 7*. The POWHEG BOX identifies these identical 
contributions initially, and stores the proportionality constants. When computing the 
Born cross section for all needed flavour structures, it computes only the minimum number 
of squared amplitudes it needs, and obtains the others using the proportionality relations 
found initially. 

2.4 The virtual amplitudes 

The user should provide a subroutine 

setvirtual (p(0 : 3, 1 :nlegborn) , vf lav(l :nlegborn) , virtual) , 

that returns in virtual the finite part V fln of the virtual cross section for the process 
with flavour structure vf lav and external momenta p. The V fln contribution is defined, in 
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conventional dimensional regularization (after renormalization) , by 



(2.11) 



L " i,j 

where the a and dj coefficients do not depend upon e. Their explicit form can be found, 
for example, in ref. [32]. The normalization factor is 



Here it is important to remark that, in the conventional dimensional regularization, B and 
Bij in formula (2.11) are evaluated in (4 — 2e) dimensions, and thus do depend upon e. In 
fact, all formulas for the soft contributions and the collinear remnants used in the POWHEG 
BOX are computed in the MS scheme, dropping the 1/e 2 and 1/e contributions written in 
terms of the (4 — 2e)-dimensional expression of the Born squared amplitudes and with 
the normalization factor of eq. (2.12) (see for example appendix A). Thus, if the virtual 
contribution is written as in formula (2.11), the divergent terms dropped in the POWHEG 
BOX cancel exactly the 1/e and 1/e 2 terms in (2.11), leaving the finite contribution V fln . 

If the NLO calculation has been performed in the Dimensional Reduction (DR) scheme, 
in order to use it within the POWHEG BOX, we do the following. First of all, we assume that 
the NLO result is expressed in terms of the standard MS coupling constant. If this is not 
the case, the appropriate straightforward corrections need to be applied. Then one defines 
VjP n R using the same formula (2.11) where now B and By are evaluated in four space-time 
dimensions. In other words, V° n R is defined as the zeroth order coefficient of the Taylor 
expansion of V° R 27r/(a s AT) in e. The expression to be used in the POWHEG BOX is [33] 



where 7(g) = N c /6, and 7(g) = (TV 2 — l)/(4iV c ), with the index i running over all coloured 
massless partons of the amplitude and where N c is the number of colours. 

The scale Q in eq. (2.12) is completely arbitrary in this context. It may be chosen 
for convenience equal to s or equal to the renormalization/factorization scale. Within the 
POWHEG BOX it has been fixed to be equal to the renormalization scale fi R , and so also the 
user should set it to this value. 

2.5 The real squared amplitudes 

Together with the Born and Born-correlated squared amplitudes, the user of the POWHEG 
BOX should provide the routine 

setreal (p(0 : 3,nlegreal) , rf lav(l :nlegreal) , amp2) 
that computes, given the momenta p of the external particles, the squared amplitude 3 for 
the real process specified by the flavours rflav, stripped off by a factor q s /(2-7t). 

3 Averaged over spin and colours; as for the Born and the virtual contributions, also in this case the flux 
factor should not be included by the user. 




(2.12) 




(2.13) 
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As for the Born contribution, also the real squared amplitude is requested for each 
individual flavour structure contributing to the real cross section. Many different flavour 
structures will return identical or proportional values for the real squared amplitude. The 
PDWHEG BOX identifies these identical contributions initially, and stores the proportionality 
constants. When computing the real cross section for all needed flavour structures, it 
calculates only the minimum number of squared amplitudes it needs, and obtains the 
others using the proportionality relations found initially. 

2.6 Analysis routines 

In the analysis file, pwhg_analysis . f , the user provides her/his own analysis routine to 
plot kinematic distributions. 

In order to have a unique analysis file, able to deal with events at the partonic or at 
the hadronic level (i.e. after the shower), we pass all the kinematics and properties of the 
particles via the HEPEVT common block. In this way, the analysis subroutine has only to 
have access to this common block and to the value of the differential cross section. This last 
subroutine is then completely independent from the program that has filled the common 
block, and can be called by PYTHIA or HERWIG too. 

The POWHEG BOX, during the integration of the B function, can produce plots with 
fixed NLO accuracy. This is done via a call to the routine analysis_driver, that fills the 
common block with the kinematic momenta. This routine receives as input parameters 
the value of the cross section at that specific kinematic phase-space point and the integer 
variable ikin. If ikin is set to 0, the event is treated as a Borndike one, and the nlegborn 
momenta in the kn_pborn kinematic common block arc copied on the HEPEVT common 
block. Otherwise, the subroutine copies on the common block the nlegreal momenta 
stored in kn_preal. As final step, it invokes the routine analysis, that receives as input 
parameter only the value of the cross section. 4 

3. The singular regions 

For each real flavour structure, the POWHEG method requires that one decomposes the real 
cross section into the sum of contributions that are divergent in one singular region only. 
In the notation of ref. [2] we write 



Each a r is thus associated with a single flavour structure, and a single singular region. 
Sometimes we will refer to it as a r region (or air region, which is the name of the variable 
that we often use to represent it in the code). The separation in eq. (3.1) is illustrated in 
detail in section 2.4 of ref. [2]. 

We have investigated two methods for obtaining such separation: the first one is based 
on the subtraction method proposed by Catani and Seymour (CS) [34], and the second 

4 In the POWHEG BOX package, we have included our own analysis routine, that uses TOPDRAWER as 
histogramming tool for our plots. 




(3.1) 
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one on the method proposed by Frixione, Kunszt and Signer (FKS) [35, 36]. It was found 
that the FKS subtraction method is better suited for our purpose, and so we focused upon 
it. The difficulties with the CS dipole method are due to the large number of dipoles and 
dipole types, and to the fact that it is not easy to separate the real cross section into a 
sum of terms having the same singularity structure of the CS dipoles. It is in fact not 
possible to simply weight the real cross sections with factors that vanish in all but one 
singular region, since for each singular region there are several CS dipoles, depending upon 
the choice of the spectator. 

The FKS method is slightly more cumbersome than the CS method when counterterms 
are computed using the collinear and soft plus-distributions. It turns out, however, that 
this difficulty remains the same for all processes. The procedure to disentangle the plus- 
distributions can be coded simply in a general way once and for all. 

In the FKS framework, singular regions are characterized as follows: 

• A final-state parton i is becoming collinear to either initial-state partons j = 1, 2, or 
soft 

• A final-state parton i is becoming collinear to a final-state parton j, or soft. 

We will call i the emitted or radiated parton, and j the emitter. If we replace the pair 
emitted-emitter with a single parton of the appropriate flavour, we obtain the flavour 
configuration of the underlying Born. 

Given the list of flavours of the real graphs, one is faced with the combinatoric problem 
of finding all singular regions. In the POWHEG framework, this task is carried out keeping 
in mind that we should be able to group easily all singular regions that share the same 
underlying Born. In order to ease this task, it is convenient to choose a standard ordering 
for the flavour structure of each air. One can easily demonstrate that the flavour of all 
the air can be ordered in such a way that the two following properties are satisfied: 

Property 1 The emitted parton is always the last parton. 

Property 2 The underlying Born configuration, obtained by removing the emitted parton, 
and replacing the emitter with a parton of the appropriate flavour, is exactly equal to one 
of the Born flavour structures present in the list of Born processes f lst_born. 

Property 2 is non-trivial, since in general the underlying Born structure obtained from a 
generic air will be equivalent up to a permutation to a flavour structure present in the 
f lst_born list. 

It is clear that putting all the air in a standard form with the properties 1 and 2 
simplifies the POWHEG implementation. In fact, since real contributions sharing the same 
underlying Born are often grouped together in POWHEG, it is better if the underlying Born 
flavour structures are unique. This will be illustrated more clearly in the example described 
in section 3.1. 

In case of initial-state collinear singularities, the emitter will be assigned the value 
1 (2) to distinguish collinear emissions from incoming line 1 (2) (the © and directions 
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of ref. [2]). If the emitted parton is a gluon, the emitter will be assigned the conventional 
value 0, that means that both 1 and 2 may have emitted it. These distinctions are such 
to minimize the number of regions, maintaining however the fact that, to each region, we 
associate a unique underlying Born configuration. 

The FKS framework, for final-state radiation (FSR) , distinguishes between the emitter 
and the emitted parton (often called the FKS parton), in the fact that only the emitted 
parton leads to soft singularities. Thus, the emitter can be a quark, with the emitted parton 
being a gluon, but not viceversa. On the other hand, the emitter-emitted can be a quark- 
antiquark pair. In this case, it does not matter what we choose to be the emitted parton. 
By convention, we will always choose it to be the antiquark. If the emitter and radiated 
partons are both gluons, when computing the corresponding a T , we supply a damping 
factor that removes the soft singularity of the emitter, with an appropriate compensating 
coefficient. For example, if we call E cm and E r the CM energies of the emitter and emitted 
parton, the damping factor may be chosen equal to E cm /(E cm + E T ), and an extra factor 
of 2 is supplied to account for the region where the role of emitter and emitted partons are 
exchanged. 

In the POWHEG BOX, the task of finding all regions associated with a given flavour 
structure is performed by the routine 

f ind_regions (nleg , rf lav , nregions , iregions) 
where nleg=nlegreal and rf lav(l :nlegreal) is the input flavour structure. It returns 
in nregions the number of regions found, and, for every found region j, it returns the 
positions (in the string of flavours) of the emitter and of the emitted parton (arbitrarily 
ordered) in iregions (1 : 2 , j ) . 

The algorithm for finding the final-state regions is the following: 

1. Loop over all massless parton pairs i=f lst_lightpart :nlegreal, j=i+l :nlegreal. 

2. Check if i and j can come from the splitting of the same parton (i.e. if they have 
opposite flavours, or if they are both gluons, or one of them is a gluon). 

3. If they cannot come from the same parton flavour, skip them. 

4. If they can come from the same parton flavour, build up a flavour list with nlegborn 
elements, obtained from rflav by deleting partons i,j and adding a parton with the 
appropriate flavour (i.e. if i , j have opposite flavour, or are both gluons, add a gluon; if 
one is a gluon add the flavour of the other parton) . Check if the newly built flavour list 
is an admissible flavour structure for the Born cross section. This is done by checking if 
the Born flavour structure at hand is equivalent, up to a permutation of the final-state 
partons, to any element of the list f lst_born. 

5. If the underlying Born flavour structure is valid, increase nregions and set 
iregions (1 , nregions) = i and iregions(2,nregions)=j. 

The initial-state regions are treated similarly. We check, for each final-state light coloured 
parton j , if it may come from the splitting of an initial-state parton. If it comes from the 
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first (second) incoming line, then the number of regions nregions is increased and we set 
iregions (1 ,nregions) = 1 (2) and iregions (1 ,nregions) = j. If the emitted parton 
j is a gluon, then only one region is generated, and we set iregions (1 , nregions) = 0. 

The first task of the POWHEG BOX is to build the list of all a T in the standard form 
specified in Properties 1 and 2. This task is performed by the subroutine genf lavreglist, 
that, more specifically, does the following: 

• It sets the variable f 1st jialr to the total number of inequivalent singular regions a r 
found. 

• It sets the variable f 1st jiregular to the number of real flavour structures that do 
not have any singular region, 5 and it fills the array 

f lst_regular (k=l :nlegreal , alr=l : f lst_nregular) 

with the corresponding flavour structures. If f lst_nregular is greater than 0, also 
the flag flg_withreg (defined in the header file pwhg_flg.h) is set to true. 

• It fills the array f lst_alr(k=l :nlegreal, alr=l : f lst_nalr) with the flavour struc- 
ture corresponding to the given air region. The ordering is guaranteed to respect 
Properties 1 and 2. 

• It sets the array f lst_emitter (air) to the emitter of the air structure. 

• It sets the f lstjmilt (air) array to the multiplicity of the alr th structure. This 
number arises because regions may be found that are equivalent by permutations of 
the external legs. Only one is retained in this case, with the correct multiplicity. 

• It sets the array f lst_uborn(k=l :nlegborn, air) to the flavour structure of the 
underlying Born of the alr th structure. 

• It sets the array f lst_alr2born(alr) to an index in the f lst_born array, that points 
to the underlying Born flavour structure of the alr th region. 

• It sets up a pointer structure from the underlying Born index to the set of air's 
that share it: it sets the array f lst_born2alr (0, jborn) to the number of air re- 
gions that have jborn as underlying Born, i.e. f lst_alr2born(alr)=jborn, and sets 
f lst_born2alr (i , jborn) , with i=l : f lst_born2alr (0, jborn) , to the correspond- 
ing air index. 

• For each air, a list of all singular regions for the corresponding flavour structure is 
also build. These are needed in order to compute 7£ Qr , as we will see further on. The 
array element f lst_allreg(l , 0,alr) is set to the number of singular regions of the 
flavour structure of the alr th region, i.e. flst_allreg(l,0,alr)=nregions of that 
particular air. Then f lst_allreg(i=l : 2 , k=l :nregions ,alr) is set to the list of 
pairs of indexes characterizing the emitter and the emitted parton for each singular 
region. 

5 An example of such configurations is given by the qq — > Hg flavour structure in Higgs boson production 
via gluon fusion. 
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In order to perform this task, the subroutine genf lavreglist loops through the f lst_real 
list of real flavour structures, calling for each of them the routine f ind_regions. Each 
region found is first transformed with a permutation, in such a way that the emitted 
parton is always the last in the list. In the case of final-state radiation, the emitter is 
also moved near the emitted parton (i.e. at the nlegreal-1 position) with a permutation. 
At this stage, one also makes sure that, if the emitter-emitted pair is made by a quark 
(antiquark) and a gluon, the gluon is always the emitted parton, and if the pair is a quark 
and antiquark, the antiquark is always the emitted parton (if this is not the case, emitter 
and emitted partons are exchanged). The lists flst_alr and f lst_emitter are updated 
accordingly. Once this procedure is completed, the air list is complete, but each element 
may appear more than once. The list is thus searched for equivalent elements, it is collapsed 
in such a way that each element appears only once, and a multiplicity factor f lst_mult 
is setup to keep track of how many occurrences of a given contribution are present. At 
the end of this procedure, only inequivalent air's remain, with an associated multiplicity 
factor. But it may still happen that the same underlying Born configuration may appear 
in different air with different ordering. At this stage the air list is scanned again. If at 
any point one finds and underlying Born that differs from one appearing in the f lst_born 
list, the flavour structure of the current arl (and its underlying Born flavour structure) 
are permuted, in such a way that one recovers the same ordering of one of the elements 
appearing in the f lst_born list. In case of final-state radiation, the emitter may end up to 
be no longer the nlegreal-1 leg of the process, and so, the array f lst_emitter is updated 
accordingly. 

3.1 Example 

Due to the intrinsic complexity of the procedure for finding the singular regions, it might 
be useful to examine it on a simple example. Let us consider a process that, at the Born 



jborn 


processes 


f lst_born 


1 

2 


sc — > gud 
gu — > ssc 


[3,4,0,2,1] 
[0,2,-3,3,4] 



Table 1: Flavour structure of two Born subprocesses and the corresponding POWHEG BOX notation, 
in the third column. 

level, has only two flavour structures (f lst_nborn=2) and five coloured massless partons 
(nlegborn=5 and f lst_lightpart=3). In the second column of table 1, we list two partonic 
Born subprocesses, with the corresponding POWHEG BOX flavour structure f lst_born. 

Suppose now that the real-process contributions are only the four ones in table 2, so 
that f lst_nreal=4. Since the real contributions have one more parton with respect to the 
Born diagrams, we have also nlegreal=6. Note that we deliberately have not included 
subprocesses such as sg — > gudc, gg — > sscu and gu — > sscg, in order to keep the example 
short. 
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jreal 


processes 


f lst_real 


1 


sc - 


> gudg 


[3,4,0,2,1,0] 


2 


sc - 




[3,4,3,-3,2,1] 


3 


gc - 




[0,4,0,2,1,-3] 


4 


du - 


->■ dssc 


[1,2,1,-3,3,4] 



Table 2: Flavour structure of four real subprocesses and the corresponding PDWHEG BOX notation, 
in the third column. 

The subroutine f ind_regions is called on each flavour structure of the list f lst_real. 
This subroutine returns the list of emitter-radiated pairs, and their total number. For 
example, when the first flavour structure [3,4,0,2,1,0] is passed to the subroutine, this 
returns the list of 7 pairs: {(3, 4), (3, 5), (3, 6), (4, 6), (5, 6), (0, 3), (0, 6)}. The (3,4) pair 
means that the gluon in the third position can be emitted by the up quark in the fourth 
position, and that, by removing this gluon, the so obtained underlying Born is a valid 
Born, since its flavour structure is equivalent (up to a permutation of the final-state lines) 
to the first flavour structure in f lst_born. The (0, 3) pair represents the singular region 
associated with gluon 3 being emitted by both initial-state partons, and, once removed, 
we get a valid flavour Born. When the subroutine is called on the second flavour structure 
in f lst_real, [3, 4, 3, —3, 2, 1], it returns a single pair (3,4), meaning that the ss pair can 
come from the splitting of a gluon, and the diagram obtained by replacing the ss pair 
by a gluon is a valid Born. Similarly the call of f ind_regions on the third element of 
f lst_real returns the pair (1, 6) that signals the fact that this diagram is compatible with 
the splitting of an initial-state gluon into an ss pair. 



air 


iregions 


f lst_emitter 


f lst_alr 


1 


(3,4) 


5 


[3,4,1,0,2,0] 


2 


(3,5) 


5 


[3,4,2,0,1,0] 


3 


(3,6) 


5 


[3,4,2,1,0,0] 


4 


(4,6) 


5 


[3,4,0,1,2,0] 


5 


(5,6) 


5 


[3,4,0,2,1,0] 


6 


(0,3) 





[3,4,2,1,0,0] 


7 


(0,6) 





[3,4,0,2,1,0] 


8 


(3,4) 


5 


[3,4,2,1,3,-3] 


9 


(1,6) 


1 


[0,4,0,2,1,-3] 


10 


(1,3) 


1 


[1,2,-3,3,4, 1] 



Table 3: List of all the 10 singular regions found up to this stage of the program, of the emitter- 
radiated pairs, of the position of the emitter, f lst_emitter, and of the corresponding flavour 
structure, f lst_alr. 

In total we have then 10 singular regions, so that, at this stage of the program, 
f lst_nalr=10. We have listed them in the second column of table 3. The flavour structure 
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of these singular regions is then saved in the array f lst_alr after a suitable permutation 
of the final-state partons, such that the emitted parton is in the last position (nlegreal), 
and, in case of final-state radiation, the emitter is in the nlegreal- 1 position. At the same 
time, the position of the emitter is recorded in the array f lst_emitter (see the third and 
fourth column of table 3). 



air 


f lst_alr 


f lst_mult 


f lst_uborn 


1 


[3,4,1,0,2,0] 


2 


[3,4,1,0,2] 


2 


[3,4,2,0,1,0] 


2 


[3,4,2,0, 1] 


3 


[3,4,2,1,0,0] 


1 


[3,4,2,1,0] 


4 


[3,4,2,1,0,0] 


2 


[3,4,2,1,0] 


5 


[3,4,2,1,3,-3] 


1 


[3,4,2,1,0] 


6 


[0,4,0,2,1,-3] 


1 


[3,4,0,2,1] 


7 


[1,2,-3,3,4, 1] 


1 


[0,2,-3,3,4] 



Table 4: List of all the 7 singular regions found, of their flavour structure, f lst_alr, of the 
multiplicity of each singular region, flstjnult, and of the underlying Born flavour, flst_uborn. 

At this stage, the list of singular regions is scanned, and, if two elements are equivalent 
(up to a permutation of the final-state partons) and have equal emitted and radiated parton, 
one of them is removed from the list and the multiplicity factor flstjnult of that singular 
region is increased. By referring to table 3, the fourth air flavour list, [3,4,0,1,2,0], is 
equivalent to the first one, [3,4, 1,0,2,0], so that it is removed and the multiplicity factor 
of the first singular region is set to 2. This is illustrated in table 4, where the fifth and the 
seventh singular regions of table 3 have been removed and the multiplicity factors of the 
second and fourth singular region in table 4 are set to 2. 

In the last column of table 4, we collect the underlying Born flavour structures, 
f lst_uborn, corresponding to each air singular region. Each of these underlying Born 
must be equivalent (up to a permutation of the final-state partons) to one of the Born in 
the list of valid Born flavour structures, f lst_born. 

At this point, we scan the list of the underlying Born flavour structures, f lst_uborn, 
and permute the flavours in such a way to obtain exactly one element of the f lst_born 
list. Correspondingly, we reorder the flavours of the corresponding air singular region. 
For example, consider the first element of flst_uborn, [3,4, 1,0,2]. This element becomes 
equal to [3,4,0,2,1] (the first element of the flst_born list), only after the exchange of 
the element in the third position with the one in the forth position (1 o 0) followed by the 
exchange of the element in the fourth position with the one in the fifth position (2 -H- 1). 
We perform the same exchanges on the first element of flst_alr, i.e. [3,4,1,0,2,0], and 
we obtain [3,4,0,2,1,0], and, after these permutations, the emitter is the parton in the 
fourth position, so that f lst_emitter=4. By performing this task on every underlying 
Born flavour structure, we obtain the second and third column of table 5. 

As final tasks, we fill the arrays with the pointers to go from an air region to its under- 
lying Born flavour structure and viceversa. The first 6 elements in the f lst_alr list have 
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air 


f lst_alr 


emi 


a2b 


nreg 


f lst_allreg 


1 


[3,4,0,2, 1, 0J 


4 


1 


7 


{(3, 4), (3, 5), (3, 6), (4, 6), (5, 6), (0, 3), (0, 6)} 


o 
z 


[3,4,0,2, 1, 0J 


5 


1 


1 


{(3, 4), (3, 5), (3, 6), (4, 6), (5, 6), (0, 3), (0, b)\ 


3 


[3,4,0,2,1,0] 


3 


1 


7 


{(3, 4), (3, 5), (3, 6), (4, 6), (5, 6), (0,3), (0,6)} 


4 


[3,4,0,2,1,0] 





1 


7 


{(3, 4), (3, 5), (3, 6), (4, 6), (5, 6), (0,3), (0,6)} 


5 


[3,4,3,2,1,-3] 


3 


1 


1 


{(3,6)} 


6 


[0,4,0,2,1,-3] 


1 


1 


1 


{(1,6)} 


7 


[1,2,-3,3,4, 1] 


1 


2 


1 


{(1,6)} 



Table 5: List of the flavours of the singular regions after they have been reordered so that the 
corresponding underlying Born is equal to one of the flavour subprocess in the f lst_born list. From 
the third column on: the position emi of the emitter (emi=f lst_emitter), the underlying Born 
flavour structure pointer a2b (a2b=f lst_alr2born), the number of singular regions nreg associated 
with that particular air region (nreg=f lst_allreg(l , , air) ) and the list of the corresponding 
emitter-radiated pairs, f lst_allreg(l : 2 , 1 :nreg, air) . 



the first Born flavour structure as underlying Born, and only the last region has the second 
Born flavour structure, so that the elements of the array f lst_alr2born are as in the fourth 
column of table 5. Viceversa, 6 air singular regions have the first Born flavour structure 



jborn 


f lst_born2alr (0 , jborn) 


f lst_born2alr 


1 


6 


{1,2,3,4,5,6} 


2 


1 


{7} 



Table 6: In the second column, the number of singular air regions that have the jborn* Born 
subprocess as underlying Born and the corresponding list of these air regions. 



as underlying Born, while the seventh has the second one, so that f lst_born2alr(0, 1)=6 
and f lst_born2alr (0,2)=1 (see table 6). The list of these singular regions is shown in the 
third column: f lst_born2alr (1 : 6 , 1)={1, 2, 3, 4, 5, 6} and f lst_born2alr (1 : 1 , 2) ={7}. 

Finally, for every flavour structure in f lst_alr, we build the list of all the associated 
singular regions. This is done by simply calling again f ind_regions on every element of 
the second column of table 5. For every air, the total number of singular regions, nreg, 
is saved in f lst_allreg(l , , air) (see the fifth column of table 5), and all the pairs of 
emitted-radiated partons is saved in the array f lst_allreg(l :2, 1 :nreg,alr) (see sixth 
column). 

4. The B function 



In order to generate an event, P0WHEG first generates a Born kinematic and flavour config- 
uration, with a probability proportional to (see eq. (4.13) in ref. [2]) 



fb 



(4.1) 
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B^(* n ) = [S(* n ) + ^(*n)] /6 + E / f^ad #(*n+l) 



where 



a r e{a r |/6} 

+ E /yO*n,e) + E /^ G e e (*n,e). (4.2) 

Q eG{ Q el/f>} aee{ a el/b} 

The square brackets with a suffix represent a context: everything inside is relative to that 
suffix. The symbol /& labels a Born flavour structure. Thus, due to this context notation, 
B and V in the square bracket refer to the contribution of the Born and soft-virtual cross 
section having the flavour structure The suffix in the first summation represents a 
sum over all a T that have /& as underlying Born. The corresponding square bracket under 
the integral sign means that we integrate in the radiation phase space of the current a T , 
keeping the underlying Born variables fixed and equal to 3> n . The real contribution R 
has been properly regularized using the plus distributions (see eq. (2.88) and (2.89) in [2]), 
so that it is written with a hat, and has to be handled properly. The way we deal with the 
generation of the underlying Born kinematics in POWHEG is the following. For each singular 
region, we parametrize the radiation variables <3? ra d as function of three variables (xrad in 
the code) in the unit cube, 



-Xrad 



'y(l) y(2) Y (3) 
^rad> ^rad' ^rad 



}■ 



(4.3) 



and we also parametrize the z variable in eq. (4.2) as a function of X^ d in the [0, 1] interval. 
We then introduce the B function, defined as 



B(* n ,X rad ) = E^ /6 (*n,*rad), 
fb 



(4.4) 



with 



x iad ) = [b(*„) + n*n)] /6 + E 

«r6{ar|/i,} 

dz 



rad 



dX 



rad 



#(*n+l) 



aes{ael/f>} 



dX, 



(i) 

rad 



+ E 7^IY G e ffi (*«,e)+ E \ 



"ee{ael/&} 



dz 



dX, 



(i) 

rad 



0* n , e ).(4.5) 



One now integrates B(& n ,X ra d) in the full (<& ra ,X ra d) phase space, using an integration 
routine that implements the possibility of generating the integrand with a uniform weight 
after a single integration. The routine that we use is mint [37], that was explicitly built 
for application in POWHEG. Since mint is designed to integrate in the unit hypercube, also 
the Born phase space has to be mapped in a hypercube, spanned by the variables X^, orn , 
as described in section 2.2. 



4.1 The btilde function 

The btilde (xx,wwwO, if irst) function implements the function B in eq. (4.4). The first 
elements of the array xx correspond to xborn and the last 3 correspond to xrad. The 
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weight wwwO is passed by the integration routine, and equals the weight factor (arising from 
integration volume and importance sampling) supplied by the integration routine. The flag 
if irst has a rather involved use (described in detail in ref. [37]) that is needed in order to 
implement the folding of some integration variable. The subroutine mint may be in fact 
requested to use more points for some selected integration variable, keeping fixed all the 
others, for each single random contribution to the integral. In practice, one may require that 
more points for the variables -X" ra( j are used for a single value of X^, orn . Upon the first call 
with a new value of -X"b orn , the function btilde is called with if irst=0. In all subsequent 
calls with the same X\> oin , but different X rax ±, btilde is called with if irst=l. After all 
calls with if irst=l, a last call with if irst=2 is performed, where btilde accumulates 
all the values computed since the last ifirst=0 call. The quantities that depend only 
upon Xbom are computed only once when the call with if irst=0 is performed. Thus, the 
Born phase space is generated at this stage. The Born cross sections for all Born flavour 
structures are computed and made available in appropriate common blocks at this stage, 
by calling the subroutine allborn. The virtual contributions to btilde are also computed 
here. The subroutines btildeborn(resborn,www) and btildevirt (resvirt , www) fill the 
arrays resborn and resvirt with the contributions for each Born flavour structure. The 
contributions from the collinear remnants and from the real cross section (btildecoll and 
btildereal) are computed both for if irst=0 and ifirst=l. Notice that, in these cases, 
the index in the arrays btildecoll and btildereal refers to a given underlying Born. 
Thus, for each array entry in btildereal, for example, several contributions from different 
a r regions (sharing the same underlying Born) are summed up. All the contributions to 
btilde are accumulated in the array results. When called with if irst=2, the behaviour 
of the program depends upon the setting of the flag negf lag. If true, btilde should only 
compute the contribution of negative weights. Thus, the positive entries in results are 
zeroed, and the negative entries are replaced with their absolute value. If negf lag is false 
(normal behaviour), the negative entries are zeroed. The array results is also stored in 
the rad_btilde_arr array, defined in the header file pwhg_rad.h. It is needed at the stage 
of event generation, where the underlying Born flavour structure will be chosen with a 
probability proportional to its entries. 

4.2 The Born cross section 

The Born contribution to the btilde function is evaluated as follows. When btilde is 
called with if irst=0, the Born phase space is generated with a call to the gen_born_phsp 
subroutine, that, in turn, calls the user provided born_phsp subroutine. For cases when 
the Born cross section is itself infrared divergent (for example, in the Z + jet production 
case), we need either a generation cut for the underlying Born configuration, or a Born 
suppression factor. The PDWHEG BOX provides a standard form for the latter. We do not 
further discuss this issue here; in the forthcoming reference [13] we will illustrate this 
problem in great detail. The factorization and renormalization scales are set with a call 
to setscalesbtilde, and then the routine allborn is called. The routine computes the 
Born contributions for all Born flavour structures, and stores them in the arrays br_born 
for the cross section, br_bornjk for the colour correlated Born cross section, and br_bmunu 
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for the spin correlated one. Many subsequent calls make use of the Born cross sections, so 
it is mandatory that this call is performed first. In particular, soft and collinear remnants 
need the Born terms, and so do the singular limits of the real cross section. 

In order to understand the code of the allborn routine, it is better to assume first 
that the flag f lg_smartsig is set to false. The role of this flag is explained in detail in 
section 4.8. 

The call to btildeborn(resborn) , in the btilde function, fills the array resborn 
with the Born cross section for each Born flavour configuration, including the corresponding 
parton distribution function (pdf) factors. In case the f lg_nlotest is set, while computing 
the integral of the btilde function, an analysis routine is also called to perform a bare 
NLO calculation for several user-provided distributions. In the case of the Born result, this 
analysis routine is called within btilde at the end of a given folding sequence. 

4.3 The soft-virtual cross section 

The soft-virtual amplitude is described in detail in section 2.4.2 of ref. [2] in the case of 
massless coloured partons. We complete here the formulae given in [2] in order to include 
the case of massive partons. All these formulae are implemented in the POWHEG BOX, that 
can also deal with massive coloured partons. When massive coloured partons are present, 
formula (2.99) in [2] becomes 



^ = ^ ( QB + J^Iij :Bij - B + V fin ] . (4.6) 



The quantity Zy is non vanishing for % and j denoting any coloured initial and final-state 
partons. X« is non zero for i denoting any massive coloured parton. They both arise from 
soft radiation, and are reported in appendix A in the same form that appears in the code. 
The Q term has the same expression given in ref. [2] 



^Eh-lo^^-C.log^) 

+2Ch ( log2 7T " log2 Cc ) " 27/1 log 7T 



- log [7/e + 2C/ e log ? c + 7/ e + 2C /e log Q . (4.7) 

Ei denotes the energy of parton i in the partonic CM frame, and fi denotes the flavour, 
i.e. g for a gluon, q for a quark and q for an antiquark. In addition we have 
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% = {% ~ 2 4- Va - 2 4t, m , i'=i'-q = [^-^r)Cp. (4.10) 



3 J 9 

We stress that now the index i in the sum of eq. (4.7) runs only over the massless coloured 
final-state partons. In fact, the contributions in square bracket arise from collinear (rather 
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than soft) final-state singularities, and thus apply only to massless partons. The last line 
arises from initial-state collinear singularities. 

The parameters 5o and £ c are arbitrary. In the framework of the POWHEG BOX, we have 
set 8q = 2 and £ c = 1. An analogous parameter for initial-state collinear singularities, Si, 
that appears in the collinear remnants, is also set to 2. These parameters are hardwired in 
the code, since there is no reason to change them. 

The soft-virtual contribution of eq. (4.6) is implemented as follows. The POWHEG BOX 
has access to the user-provided Born cross section, including colour correlations. It thus 
builds automatically the Q, lij and Xj contributions of eq. (4.6) to the soft-virtual cross sec- 
tion. The corresponding code is found in the btildevirt subroutine, in the sigsof tvirt . f 
file. As already stated in section 2.4, the scale Q is chosen equal to the renormalization 
scale fi R . The user-provided virtual cross section should then have Q = fj, R . The subroutine 
btildevirt fills the array resvirt with the contribution of the soft- virtual cross section 
for all possible underlying Born configurations. 

4.4 The collinear remnants 

The collinear remnants G^® and Gq 6 of eq. (4.2) are the finite leftover from the subtraction 
of collinear singularities. Their general form, in the FKS framework, is given in eq. (2.102) 
of ref. [2]. They are implemented in the POWHEG BOX in the MS scheme. 6 Their implemen- 
tation in the btilde function has to properly handle the distributions in the z variable. 
This is done using the identities 

1 ' +iogi^«y(i-z), (4.11) 
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where x stands for either x<$ or xq. The z integration extends from x to 1, and thus it is 
performed using the rules 

\z f{z)-^— = fdz + log /(l), (4.13) 

1 /l — ,(-1 ~\\ /•! 1 — ,fi ~\ 1 — ,2/i \ 1 — 2 



l dz / w (ML^£)) _jf i, _ /(1 ),M1Z£) + ^d-,)-^6, /(1)| (4 , 4) 



x \ ± * /(, 

where we have always set £ c = 1 in the code. In this case too, a loop over all underlying 
Born configurations is performed, and, for each of them, all singular remnant contribu- 
tions appropriate to the flavours of the corresponding incoming legs are computed. The 
btildecoll function requires an extra integration variable, given as its first argument, 
to generate a value for z. Within the btildecoll function, if the f lgjilotest flag is 
set, a call to the analysis routines is performed to output the contribution of the collinear 
remnants to the user-defined kinematic plots. 



"Since currently the DIS factorization schemes are no longer used, we do not implement them. 
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4.5 The real contribution to btilde 



The real contribution to the btilde function has a certain complexity, mostly due to the 
handling of the distributions in the £ and y variables. The function btilde simply calls the 
function btildereal to get the contribution of the real cross section for each underlying 
Born configuration. The complex task of evaluating the real integrand is carried out in 
the subroutine btildereal. In this subroutine, there is a loop over all possible emitters, 
that envelopes its whole body. All contributions are then evaluated for a given emitter. 
With this choice we avoid repeating continuously complex phase-space evaluations. For 
each emitter, for given underlying Born and radiation variables, there is only one real 
kinematics to consider. 

If the emitter is a final-state parton, the real phase space is generated by a call to the 
subroutine gen_real_phsp_f sr. For initial-state emission, a call to gen_real_phsp_isr is 
done. The program then calls the subroutine sigreal_btl, that fills its array argument 
with the contributions of all regions (i.e. air) that have as emitter the current emitter 
kn_emitter. This subroutine returns an array whose elements are 

in the FSR case and i = 2 in the ISR case, rather than R a alone (R a is defined in the 
notation of [2]). The quantity R a (l — y 1 )^ 2 has well defined soft, collinear and soft-collinear 
limits, that are obtained by a call to the subroutines soft, collf sr (for final-state radia- 
tion) or collisrp and collisrm (for initial-state © and © collinear regions), sof tcollf sr 
for soft-collinear limit in final-state radiation, and softcollisrp and sof tcollisrm for 
soft-collinear initial-state radiation. The handling of the £ and y distributions require some 
care, and we describe it here in some detail. 

We begin by looking at the final-state radiation case. The integral that we would like 
to perform has the form 



the Jacobian being given by formula (5.40) in ref. [2]. We have indicated explicitly the 
dependence of the Jacobian upon the radiation variables, but one should keep in mind that 
it also depends upon the underlying Born variables. We have assumed, in all generality, 
that the upper limit in the £ integration in eq. (4.15) may depend upon y, and we have 
denoted it as X(y). This is in fact not the case in our choice of the final-state radiation 
kinematics, but it is the case for initial-state radiation (ISR). The gen_real_phsp_f sr 
subroutine returns the Jacobian for the integration of the radiation variables divided by £, 
in the variable jac_over_csi. 

In eq. (4.15) we introduce a new rescaled variable £, whose upper bound does not 
depend upon y 




[(i - y) ^ ] 6) + (i^) + ' (4 - l5) 



(4.16) 



We can easily show that 




(4.17) 
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We now rewrite eq. (4.15) as 
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where we should not forget that £ is now also a function of y through eq. (4.16). Defining 



f(Z,v) 



[(1 - 2/K 2 i? Q ] , 



(4.19) 



we get 



+ [lo g X(y)/(0,y)-logX(l)/(0,l)] 



f(SX(y),y) - f(0, y) /(f 1) - /(0, 1) 



(4.20) 



We now see that both (1 — y) ^ 2 R a and J(£, y, </>)/£ in eq. (4.19) should be computed also 
with £ = (the soft limit), y = 1 at fixed £ (the collinear limit) and both £ = and y = 1 
(the soft-collinear limit). 

The case of initial-state radiation is handled similarly, except that now our starting 
formula is 
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(4.21) 



and one proceeds as before, by treating the y = 1 and y = —1 regions independently. 
4.6 The btildereal subroutine 

We now give a more detailed description of the way btildereal is implemented. First of 
all, the generation of the radiation phase space is performed according to the description 
given in section 5.1.1 of ref. [2] for initial-state radiation, and in section 5.2.1 for final-state 
radiation. The subroutines gen_real_phsp_f sr and gen_real_phsp_isr generate the phase 
space as a function of the underlying Born kinematics, and of three real variables xrad(3) , 
that assume random values between zero and one. Since we use eq. (4.20), a common 
Jacobian factor xjac for the transformation xrad(3) into £, y and (j) is also provided. The 
program also sets the variables jac_over_csi, jac_over_csi_coll and jac_over_csi_soft 
to J(£, y, 4>)/£, to its collinear limit and to its soft limit respectively, and multiplies them 
by xjac. 

In btildereal, these limits are obtained through the calls to soft, collfsr and 
softcollfsr, and the limits for J(£, y, (/))/£ are provided by the phase space subroutines 
gen_real_phsp_f sr, in the variables jac_over_csi_coll and jac_over_csi_soft. Notice 
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also that, while the soft limit of </(£,y, is y independent, this may not be the case 
for (1 — y)^ 2 R a ■ The computed values of the real contribution, the soft, collinear and 
soft-collinear counterterms are divided by (1 — y) £ (as in eq. (4.20)) and accumulated with 
the appropriate sign in the array resreal, indexed by the underlying Born index of the 
current air. In the case of final-state radiation, the upper limit for £ does not depend 
upon y, being given by formula (5.49) of ref. [2], and is set by the phase space program 
gen_real_phsp_f sr in the common variable kn_csimax. Its value is taken from the array 
kn_csimax_arr, indexed by kn_emitter, that is filled by the routine gen_born_phsp when 
the underlying Born phase space is generated in btilde. In case of initial-state radiation, 
kn_csimax is y dependent, and is computed by an appropriate routine when the real phase 
space is generated. The last two terms in eq. (4.20), ^-independent, are also accumulated 
in the resreal array. 

After the appropriate calls to the routine that generates the real contributions and 
its various limits, the program loops through all real regions (i.e. air), and accumulates 
the real contributions according to eq. (4.20) or its initial-state radiation version, in an 
array indexed by the index of the underlying Born of the current air. This was set in the 
combinatoric routines, in the array f lst_alr2born. The accumulated values include the 
underlying Born Jacobian, the real contribution or one of its various limits, jac_over_csi 
or one of its limits, and the remaining factor of 1/(£(1 — y)) for FSR, or 1/(2£(1 ± y)) in 
the ISR case. The sign of each contribution can be read out from eq. (4.20). Besides filling 
the output array resreal, if the flag f lgjilotest is set, the result is also output to NLO 
analysis routines, that perform a parton-level NLO calculation of user-defined distributions. 
The analysis driver for the NLO output, i.e. the subroutine analysis_driver, described 
in section 2.6, is then called with the flag set to 1 for the true real contribution, and for 
all remaining terms. 

4.7 The subroutine sigreal_btl 

The subroutine sigreal_btl fills its output array argument, indexed by the air, with the 
real contributions that have as emitter kn_emitter. The real contribution should also be 
multiplied by £ 2 (1 — y) for FSR, or £ 2 (1 — y 2 ) for ISR, and, furthermore, should also be 
multiplied by the S a functions, described in section 2.4 of ref. [2]. Two flags control the 
behaviour of this function: f lg_smartsig and f lg_withdamp. An explanation of the role 
of f lg_smartsig is given in section 4.8. The f lg_withdamp flag will instead be explained 
in section 5. 

The code is better understood if one assumes, to begin with, that these flags are set 
to false. In this case, the program loops over all air. For those that have emitter equal 
to the current emitter, it calls the subroutine realgr, passing as argument the list of 
flavours of the current configuration, and the real momenta. The function is supposed to 
return, in its last argument, the value of R, the real matrix element squared. Next, each 
contribution should be multiplied by its S a factor. In the framework of the POWHEG BOX, 
we define the S a factor in the following way. We consider the flavour structure of the a 
region under consideration, and call it f a . We call Rf<* the corresponding real contribution. 
Rja can have several singular regions (see, for example, the list of pairs of indexes in the 



- 25 - 



last column of table 5). Each such region is characterized by a pair of indexes in the legs 
of the real process, of the form These can be the indexes of two final-state lines 

becoming collinear, or of an initial- and final-state line becoming collinear. According to 
the POWHEG conventions, one can also set % = 0, meaning that there are initial-state collinear 
singularities in both directions (gluon emission from initial-state partons), and they share 
the same underlying Born. We call X a the array of all singular regions, i.e. X a is an array 
of pairs of indexes. In particular, the emitter associated with the a region, together with 
the last parton (i.e. the nlegreal parton) form a pair that belongs to X a . Let us call it 
the (k, n) pair. Then, S a is given by 



S a = -r- "H - (4-22) 




where dij are appropriate kinematic functions that vanish when lines i and j become 
collinear. The choice of the dij function implemented in the POWHEG BOX can be found in 
the compdij subroutine, in the gen_real_phsp . f file. When the phase space is generated, 
compdij is called, and the array kn_dijterm is filled. It is an ordered array, i.e. one always 
assumes i < j. For initial-state singularities it is given by the expressions 

d 0j = [E](l-y])] p \ (4.23) 
d,j = [2E](l- yj )] Pl , (4.24) 
d 2j = [2E](l + yj )] pi , (4.25) 

where yj is the cosine of the emission direction of parton j, in the real CM frame, relative 
to the positive collision direction. Notice that we assume, by convention, that k = means 
that there are collinear singularities from both the positive and negative direction with the 
same underlying Born configuration, as is the case when a gluon is emitted. The positive 
and negative collinear directions need to be considered separately if the corresponding un- 
derlying Born differs. The parameter p\ corresponds to the parameter par_diexp, defined 
in the pwhgJm.h include file. Its default value is 1. For FSR regions we define 



dij — 
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(Ei + Ej) 2 



P2 

(4.26) 



where p 2 corresponds to par_dijexp, and is set to 1 by default. 

The sigreal_btl subroutine builds the S a factor for each non vanishing air. At the 
combinatoric stage, for each air, a list of the singular regions associated with its flavour 
structure was built. This list was stored in the f lst_allreg array, and is used to find the 
indexes ij for all the singular regions of the given air (see table 5 for an example). 

One extra factor is supplied for final-state singularities if both the emitter and radiated 
partons are gluons. One multiplies the result by 



2-E'cm 
E C m + E r ' 



(4.27) 
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where E em and E r are the energy of the emitter and of the radiated parton, evaluated in 
the partonic CM. This does not change the cross section, because of the symmetry in the 
exchange of the two gluons, but guarantees that only when the radiated gluon becomes 
soft we can have a soft singularity As a final step, the multiplicity of the current air and 
the £ 2 (1 - y 2 ) (for ISR) or £ 2 (1 - y) (for FSR) factors are included. 

4.8 The f lg_smartsig flag 

In several processes, organizing the program in the way described in the previous section 
would lead to several calls to the same matrix elements, with a consequent waste of com- 
puting time. In the process of W production, for example, the matrix element is the same 
whether it is a ud or a cs collision. Or it may differ only by a Cabibbo-Kobayashi-Maskawa 
matrix element factor. If the flag f lg_smartsig is set to true, upon the first call to the 
sigreal_btl subroutine, the routine finds matrix elements that differ only by a constant 
factor, and builds appropriate array of pointers and proportionality constants. Upon sub- 
sequent entries in the program, multiple calls to proportional matrix elements will thus 
be avoided. All subprograms that invoke user routines for matrix element calculations are 
affected by the setting of f lg_smartsig, and implement the same mechanism for avoiding 
useless calls to user routines. By closely examining the code, the reader can find out how 
this works. The output of the program, however, should be independent on the setting 
of this flag. Only speed will be affected. In fact, the random number generator, used to 
set up the random kinematics to check matrix elements for proportionality, is reset to its 
original value after the equivalent matrix elements are found. 

In order to understand how the programs operate when the f lg_smartsig flag is 
turned on (i.e. is set to true) it is better to examine the allborn routine. Upon the first 
call to allborn, the current random number is saved, and then the Born cross section 
for all flavour components is computed, for several value of randomly chosen external 
momenta. An integer array equivto is set up, its value being -1 by default. If the Born 
contribution for the j th Born flavour configuration is found to be proportional to a previous 
k th Born flavour configuration (with k < j), the value of equivto(j) is set equal to k, 
and the array of real numbers equivcoeff (j) is set to the proportionality constant. Upon 
subsequent calls to allborn, this information is used to avoid further calls to setborn, 
whenever possible. Notice the use of randomsave and randomre store. By enclosing a set 
of instructions between a randomsave and a randomrestore call, we make sure that the 
random number sequence is not altered by the inserted instructions. 

If the flag f lg_smartsig is set to false, all calls to the matrix element routines are per- 
formed, but, thanks to the saving and restoring of the random number sequence, the out- 
put of the program should be independent of it. In other words, using f lg_smartsig=true 
should only accelerate the program, without altering the output. This feature can be used 
to check that nothing weird has happened in the setup phase of the equivto and equivcoef 
arrays. 

4.9 The soft, collinear and soft-collinear limit functions 

These limit functions could be obtained, in principle, by numerical methods, using the full 
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real contribution. We have, however, preferred to compute them using the factorization 
formulae for collinear singularities, and the eikonal formulae for soft emission, to avoid 
numerical instabilities. Furthermore, the real contributions, and all the manipulations 
performed by the combinatoric package can be tested for consistency (see appendix F.l). 

These soft, collinear and soft-collinear routines are collected in the file sigcollsof t . f . 
The subroutines relevant for the computation of the btilde function are reported in table 7. 
The basic formulae for the collinear limits are collected in appendix B. Here we illustrate 



soft 


soft limit 


collf sr 


collinear limit for FSR 


sof tcollf sr 


soft-collinear limit for FSR 


collisrp 


collinear in the © direction 


sof tcollisrp 


soft-collinear limit in the © direction 


collisrm 


collinear in the © direction 


sof tcollisrm 


soft-collinear limit in the © direction 



Table 7: Subroutines for the soft and collinear limits of the real contributions used in the compu- 
tation of btilde. 



the code of the collf sr routine. This routine in turns calls collf srnopdf, and provides the 
luminosity factor to its output. Another ingredient that is necessary to build the collinear 
limit functions is the direction of the transverse momentum k T of the radiated parton with 
respect to the emitter in the collinear limit (see appendix B), defined in the partonic CM 
frame. This is a function of the emitter direction and of the azimuthal angle (ft. The origin 
of the azimuth, i.e. the plane along which (ft = (or n) should be consistent with what 
gen_real_phsp_f sr does in the collinear limit. A change of (ft — > (ft + ir is instead irrelevant. 
Our convention for the origin of (ft is to take a plane containing k era (the momentum, in 
the CM frame, of the parton that will undergo the splitting in the underlying Born, see 
ref. [2]), and the third axis. The subroutine buildkperp, called from the collf srnopdf 
routine, constructs the 4- vector kperp(0:3). This vector is normalized arbitrarily, and has 
zero time component. Its only requirement is that it should be parallel to k T . Its modulus 
squared is also returned by the buildkperp routine. For each air sharing the current 
emitter, the subroutine collf sralr is called, with kperp also passed as an argument. The 
values of £ and x, defined as 

k° 2k° lt , 

and x/£ are all passed to the subroutine, that is meant to work also if £ and x vanish, with 
their ratio remaining finite. The subroutine collf sralr implements the formulae given in 
appendix B in a straightforward way, multiplying them by £ 2 (1 — y), and taking care to 
use the x/£ variable when necessary, in such a way that one never divides by x or £. There 
is only one caveat to keep in mind. The collinear approximation is meant to reproduce an 
R a contribution. In the collinear limit of the a region, this coincides with R, since the S a 
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factor becomes 1 in this limit. We should remember, however, that, in case the emitter 
and radiated partons are both gluons, we have also supplied a factor 2E ern /(E em + E T ), 
which becomes 2(1 — x) in the collinear limit (see eq. (4.27)). In addition, in this case, we 
supply a factor of 1/2 to account for the two identical partons in the final state. Thus, an 
extra factor of (1 — x) is supplied. 

The case of initial-state collinear singularities is handled similarly. In this case we 
simply have x = 1 — £, and, as before, one should evaluate all contributions taking care of 
never dividing by 1 — x. The routine collisralr implements the formulae in appendix B. 
It carries an integer argument i that corresponds to 1 for the collinear © direction and 2 
for the © direction. 

In order to get the soft-collinear limits, the subroutines softcollf sr, sof tcollisrp 
and sof tcollisrm simply call the corresponding collinear subroutines setting temporarily 
kn_csi equal to zero. 

The soft limit is obtained with the subroutine soft, that in turn calls softalr to get 
the soft contribution of a single air. It implements formula (A.l) (for e = 0), multiplied 
by £ 2 (1 — y) for an FSR region, or £ 2 (1 — y 2 ) for an ISR region. In formula (A.l) there are 
two powers of the soft momentum k in the denominator. We then factorize k° in front of 
the four-vector k and define k = k°k, so that (£ = 2fc / v / s) 




(4.29) 



is finite. The vector k carries the information of the direction of the radiated parton. In 
practice, we thus replace A; — > k in eq. (A.l), and supply the factor 4(1 — y)/s or 4(1 — y 2 )/s. 
The vector k should equal k/k° in the limit £ — > 0, keeping y and 4> fixed, for the given 
underlying Born kinematics. It is computed in the subroutine gen_real_phsp_f sr and 
gen_real_phsp_isr, with a call to the subroutine setsoftvec. Furthermore, we should 
remember that the functions S a have non-trivial soft limits. They should be computed in 
the soft limit and multiplied by the result. For this purpose, the routine compdijsoft, 
in the gen_real_phsp . f file, computes the soft limit of the dij functions, and stores them 
in the array kn_dijterm_sof t. This array has a single index, since the second one is the 
index of the soft parton. There is no need to consider the other dy terms (those not 
involving the soft parton) since in the soft limit they are finite, and do not contribute to 
S a . Observe that compdijsoft assumes that the dij terms are homogeneous in k , which 
is the case if the two parameters par_diexp and par_di j exp are equal. If they differ, the 
fastest vanishing one (in the k° — > limit) will dominate S a . It is in principle possible 
to experiment with settings that have different par_diexp and par_dijexp, provided that 
the slowest vanishing ones are excluded in some way from the list. At the moment, this 
possibility has not been investigated. As a last point, special treatment was required for 
the FSR collinear limit of two outgoing gluons. Here, in fact, no action should be taken: 
one should provide a factor (1 — x), that equals one in the soft limit. 
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5. Tuning the real cross section in POWHEG 

In POWHEG it is possible to tune the contribution to the real cross section that is treated 
with the Monte Carlo shower technique. This was pointed out first in ref. [1], where the 
POWHEG method was formulated, and it was first implemented in ref. [10]. In POWHEG there 
is the possibility to separate the real cross section, in a given singular region a, as follows 

R a = R a s +R a fj (5.1) 

where RJ has no singularities and only R" is singular in the corresponding region. In 
practice, the separation may be achieved, for example, using a function of the transverse 
momentum of the radiation ^ F{k^) ^ 1, that approaches 1 when its argument vanishes, 
and define 

K = R a F{k 2 T ) , (5.2) 
RJ = R a [l-F(k^)]. (5.3) 

One carries out the whole POWHEG-style generation using Rf rather than R a . The contribu- 
tion R^, being finite, is generated with standard techniques, and fed into a shower Monte 
Carlo as is. 

More generally F can be chosen as a general function of the kinematic variables, 
provided it approaches 1 in the singular region. This turns out to be useful in all cases 
when the ratio R/B in the POWHEG Sudakov exponent becomes too large with respect to 
its corresponding collinear or soft approximation (see for example ref. [7]). In this case, 
radiation generation becomes highly inefficient. A general solution to this problem (which 
has already been implemented in refs. [14] and [13]) is to chose the function F in the 
following way: if the real squared amplitude (no parton distribution functions included), in 
a particular singular region, is greater than five times its soft and collinear approximation, 
then F is set to zero, otherwise is set to one. We also stress that this procedure remedies 
automatically to the Born zeros problem examined in ref. [7]. 

This feature is implemented in the POWHEG BOX. By setting the flag f lg_withdamp to 
true, this behaviour is turned on. When computing the btilde function, the real contribu- 
tion will always be multiplied by a damping factor, supplied by the routine bornzerodamp. 
The damping factor is not necessary in the soft and collinear counterterm contributions, 
since, in these cases, we will certainly have F = 1. The routine bornzerodamp takes as 
argument the a-region index (i.e. the air), the value of 1Z a (i.e. the real cross section 
without the parton distribution function factors) and the value of its collinear and soft 
limits (also without pdf factors). It returns the damping factor as its last argument. The 
presence of the collinear and soft limits of lZ a in the arguments of the subroutine, allows 
the user to set a damping factor that depends upon the distance of the real contribution 
from its collinear or soft approximation, as stated previously. This routine can be easily 
modified by the user: for example, the sharp theta function adopted in the POWHEG BOX 
can be replaced by a smoother function, the factor of five can be changed, and so on. 

One can see the effects of setting the f lg_withdamp flag in the sigreal_btl subrou- 
tine. The soft and collinear limits of the real contribution are obtained with calls to the 
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subroutines collbtl and sof tbtl (the btl ending standing for btilde), that make use of 
the subroutines already described previously for the calculation of the soft and collinear 
limits. 

If a damping factor is used, the leftover term RJ of eq. (5.3) needs to be handled 
independently. The subroutine sigremnants deals with this term together with the real 
terms that do not have any associated singular region, if there are any. It has the same 
calling sequence of btilde. It is meant to be integrated using the mint integration program, 
that allows for the possibility of generating phase space kinematics distributed with a 
probability proportional to the integrand, after a single integration. Within sigremnant, 
the contribution from the regular real graphs can be integrated with an arbitrary phase- 
space parametrization, that we choose to be the initial-state radiation parametrization, 
i.e. the gen_real_phsp_isr subroutine. Both the underlying Born configuration and the 
real phase space are generated within sigremnant. The regular contributions to the real 
cross section are returned by the subroutine sigreal_reg. Within sigreal_reg, by making 
use of the list of regular real contributions (that is setup when the combinatorics is carried 
out), the regular contributions to the cross section are computed. The contribution of 
the RJ terms is more delicate. This is computed with a loop through all possible emitter 
values using the global variable kn_emitter. The real phase space is set according to it. 
Then the subroutine sigreal_damp_rem (where dampjrem stands for damp remnants) is 
invoked. This subroutine is very similar to the sigreal_btl subroutine. For all air that 
share the current emitter, the corresponding R a is computed, the damping factor dampf ac 
is computed, and the real result is multiplied by (1-dampfac) (in sigreal_btl the result 
was instead multiplied by dampf ac). 

Notice that sigreal_damp_rem and sigreal_btl carry out very similar tasks, the only 
difference being the presence of the factor (1 — F) in the first, and F in the second. 
This fact is exploited in the POWHEG BOX by implementing both of them via a call to a 
single subroutine sigreal_btlO, that carries an extra integer argument. When the extra 
argument is zero, the multiplication factor is set equal to F, and when it is 1, it is set to 
(1-F). 

6. The initialization phase 

The preparation of the grids for the generation of the events is carried out in the subroutine 
bbinit. Its most important task is to execute the integration of the btilde function, 
determine the fraction of negative weights, compute the total cross section, and, if required, 
plot the NLO distributions. At the first step, the subroutine mint [37] is invoked with 
imode set to 0. In this mode, mint integrates the absolute value of btilde, and sets up the 
importance-sampling grid. Next, mint is invoked with imode set to 1, and the flag negf lag 
set to true. In this mode, mint computes the negative contribution to the btilde function. 
No histograms for the NLO results are generated up to now. At this stage, negf lag is set 
to false, the f lg_nlotest is set to true, and mint is invoked again on btilde to compute 
the positive contribution to the integral. At this stage, the NLO histograms are filled. 
We stress that also negative weights, if present, will end up in the histograms, so that the 
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NLO histograms should exactly correspond to a standard NLO calculation. The positive 
weight total cross section computed by mint is combined with the negative weight part, 
and stored in a variable racLsigbtl, defined in the header file pwhg_rad . h. After this, the 
contribution to the cross section from the real remnants is also computed. These are terms 
that arise either because there are real contributions with no associated singular regions, or 
because f lg_withdamp is set to true (see section 5). The remnant cross section calculation 
is performed with an independent set of grids. Also the remnant contributions will end up 
in the NLO histograms. The remnant cross section is stored in the variable racLsigrm, 
and the total in rad_sigtot=rad_sigrm+rad_sigbtl. 

When mint is called with imode equal to 1, the upper bounding envelope of the inte- 
grated function is also computed, and stored in an array. This upper bounding envelope 
will be used later for the generation of unweighted events. The arrays xgrid, ymax, xmmm 
are all necessary for the generation of the events, and they can be saved in a file, so that the 
time consuming initialization phase does not need to be repeated if one wishes to generate 
more events in the same context. 

The final important task of the bbinit routine is the call to the do_maxrat subroutine, 
that sets up the normalization of the upper bounding function for radiation, thus prepar- 
ing the system for the generation of full events. This will be described in section 7.1.1. 
In bbinit an initialization call to the function gen, that generates the underlying Born 
configuration, is also performed. 



7. The generation of radiation 



There are two components that contribute to the generation of radiation: one arises from 
the B term and the other from the remnant. The total cross section for the two contri- 
butions is stored in the global variables rad_sigbtl and rad_sigrm. When radiation is 
generated, one begins by picking one of the two cases with a probability proportional to 
the respective cross section. In the PDWHEG BOX, the generation of radiation is carried out 
in the subroutine pwhgevent, that begins precisely by performing this random choice. We 
describe, in turn, the two components. 

7.1 Radiation from the B component 

This begins with the generation of an underlying Born configuration distributed according 
to the B function. Radiation is generated using the POWHEG Sudakov form factor (see 
eq. (4.21) of ref. [2]) 

A /l (*»,Pr)= II A&(*n,j>r), (7.1) 
a r e{a r |/b} 

where 



A^(* n ,p r ) = expJ- J 



- n — 



(7.2) 



If R has been separated into a regular and singular part, according to eq. (5.1), only the 
singular part will appear in the Sudakov form factor. According to the notation of ref. [2], 
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the square bracket with the a T suffix indicates that all quantities inside the bracket should 
be taken relative to the a r region. So, the n + 1 phase space in eq. (7.2) is given as 
a function of the underlying Born phase space #° r , taken at the point <fr n , and of the 
radiation variables <fr ra d, according to the phase-space mapping defined for the a T region. 
In POWHEG, the individual Sudakov form factors for each a T are further assembled into a 
product of form factors sharing the same underlying Born and the same radiation region. 
As we have seen, each singular region is characterized by an emitter and an emitted parton. 
Within POWHEG, the emitted parton is always the last one, while the emitter can be any 
light coloured parton in the initial or final state. There is one single initial-state radiation 
kinematics, independent of which incoming parton is emitting. The final-state radiation 
kinematics depends instead upon the index of the emitter. We introduce here the label 
rr to specify the radiation region kinematics: rr = 1, if kn_emitter=0, 1 or 2, and rr = 
kn_emitter - f lst_lightparton+ 2, if kn_emitter > f lst_lightparton. In fact, within 
the POWHEG BOX framework, the radiation kinematics is the same for kn_emitter=0, 1 or 
2. We write 

A^(*n,P T )= [J A rr(*n,P T ), (7.3) 
rre{rr|/ 6 } 

where 



A&(*n,J>r) = expi- J d$ tai ^Sgr 0(k T (* n+l ) - p T ) 



a r e{a r |/ 6 ,rr} 



-n — 



, (7-4) 



and the notation {cx T \fb, rr} indicates the ensemble of all a T that share the same underlying 
Born fb and the same radiation region kinematics rr. It makes then sense to define 

R"(* n+1 )= Yl Rar (*n + i), (7.5) 

a r e{a r \f b ,TT} 

since the phase space only depends upon the radiation region kinematics rr, and not on 
the specific a r . With this definition we have 



A£(*„,pr) = exp<^ - 



fi >ri Y<f> 



- n — ^ n 



(7.6) 



In order to generate the radiation, the POWHEG BOX uses the highest-bid algorithm. For 
each rr, it generates a p T value with a probability distribution equal to 

p/Hpt) = ^-A£(* n>Pr ). (7.7) 

The program then selects the highest p T value, and thus fixes the corresponding rr region. 
The a T value is picked in the ensemble {a r |/b,rr}, with a probability proportional to the 
corresponding R ai . 

The individual p T values for each Arr are generated with the veto method. We define 

d$ rad = J^dtdydfr (7.8) 
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where J rr is the Jacobian of the ^rad phase space, when written as function the three 
radiation variables: £, y and <p. We then introduce a suitable upper bounding function 
U"{i,y), and determine its normalization N" by requiring 

for all 3> n and $ ra d- In order to generate the radiation, one first generates a p T value 
according to the probability distribution 

P%(Pt) = ^K(Pt), (7.10) 



where 



A tt(Pt) = exp 



-Nf b J d^dyd^U"(ty)0(k T (^ n+1 )-p T ) 



(7.11) 



and then generates the corresponding values for the radiation variables £, y and <p, dis- 
tributed with a probability proportional to U TT . At this point one accepts the event with 
a probability 

1 J"R"(* n+1 ) 
Nf b U"(^y) Bf»(* n ) ■ [ ■ ) 

If the event is rejected, one goes back to the beginning, and generates a new p' T value, 
smaller than the current one, using the probability 

P "^ = W,W>' A(p ^ (7 ' 13) 

Ideally, the function U" should be chosen in such a way that the equation r = A^(p T ) can 
be easily solved for p T , and a set of radiation variables, distributed according to U rT , are 
easily generated. In practice we only require this second feature, and solve for the equation 
r = A^(p T ) numerically (in this way, if r is a uniform random number between and 1, 
the corresponding p T is distributed according to eq. (7.10)). Several choices are possible. 
In appendixes C and D we describe the functions used in the POWHEG BOX. 

In the POWHEG BOX, a straightforward variant of the veto method is used several times. 
Suppose that we know a function F($ Ta d) such that 

jit f>??((h \ 

Then we can first use the veto method accepting events with a probability 

-^(^rad) (7 1<\ 

Nf b u^ rad y { • ' 

If the event is accepted, we go through a second veto, accepting the event with a probability 

1 gg^+i) (716) 
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This is obviously the same as accepting according to the probability of eq. (7.12), but it 
has the advantage that, in many instances, when a veto is imposed, one only evaluates the 
function F, without the need to compute the real and Born contributions. This method is 
also applied in the PDWHEG BOX by replacing eq. (7.9) with 

TIT f>?T(<ff \ 

Bfb(&n) ^ ^ V) U " ( ^ d) ^ N Z ^ r($rad) ' (7 - 17) 

where iVj£(£,y) is a stepwise function in the £ and y radiation variables, and where 

N}l = m a xN}l(Z,y). (7.18) 

One determines the step function 7Vj£ (£,£/) so to have the smallest values that satisfy the 
first bound of eq. (7.17). Then, the method described above is used, where, according to 
eq. (7.14), 

F(<S> iad ) = N}l(Z, y )U"(^ d ), (7.19) 

so that events are first accepted with a probability Nj?(£,y)/Nj?. 

Within the POWHEG BOX, the generation of the underlying Born kinematics is per- 
formed by the routines gen_btilde, that invokes gen with the appropriate arguments. 
After that, the subroutine gen_uborn_idx is called and it generates the underlying Born 
flavour configuration. The purpose of this call is to pick a random /& configuration, with 
a probability proportional to its contribution to the B value at the given kinematic point. 
By inspecting the gen_uborn_idx subroutine, we see how this task is performed. We recall 
that when gen returns, the last call to the btilde function has been performed at the 
generated Born kinematics configuration. The contribution of each flavour component of 
the B cross section is stored in the array rad_btilde_arr. In gen_uborn_idx a generic 
utility subroutine pick_random is invoked with arguments f lst_nborn, rad_btilde_arr 
and rad_ubornidx. The pick_random subroutine returns the value rad_ubornidx with a 
probability proportional to rad_btilde_arr (rad_ubornidx) . The variable rad_ubornidx 
represents the index of the currently generated underlying Born configuration. There are 
several other rad_ prefixed global variables that need to be set, in order to perform the 
generation of radiation from the current underlying Born. First of all, a list of all air, 
that share the current underlying Born structure, should be constructed. This is done by 
filling the array rad_alr_list, of length rad_alr_nlist, using f lst_born2alr, that was 
constructed at the combinatoric stage of the program. 

The variable denoting the singular region rr is represented by the global variable 
rad_kinreg in the POWHEG BOX. As already stated, it takes the value 1, for initial-state 
radiation, and the value rad_kinreg = kn_emitter-flst_lightparton+2 for final-state 
radiation. Not all values of rad_kinreg may be associated with an active radiation region 
for the given underlying Born. A logical array rad_kinreg_on is set up, with its entries 
indexed by the rad_kinreg values. The entries set to true correspond to active radiation 
regions. The array rad_kinreg_on is set in the subroutine gen_uborn_idx. In table 8 we 
summarize the global variables relevant to the generation of radiation. We do need all these 
variables, because we typically need to consider the air that share the current underlying 
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racLubornidx 


index in the current underlying Born flavour structure 


rad_alr_list 


list of air's that share the current underlying Born 


rad_kinreg_on 


marks the active singular regions for the current underlying Born 


radjcinreg 


current singular region 



Table 8: Global variables that characterize the generation of radiation for the given underlying 
Born configuration. 



Born and the current kinematic region. This is done by going through the rad_alr_list, 
and considering only the air's whose emitter is compatible with the current radJcinreg 
value. 

7.1.1 Normalization of the upper bounding function 

Before the radiation is generated, the normalization of the upper bounding functions should 
be computed. This task is carried out by the subroutine dojnaxrat, which, in turn, is 
invoked in the bbinit subroutine. The normalizations iVj£(£,y) and Nj* are stored in the 
arrays 

rad_csiynorms (radjicsinorms , rad jiynorms , radjikinreg, flstjiborn), 
rad_norms (rad_nkinreg, flst_nborn). 

By inspecting the do_maxrat routine, one can see that there is a mechanism (that is better 
understood by studying the code) to store and retrieve previously computed values for 
these arrays. The core of the do_maxrat routine is a loop repeated nubound times, where 
nubound is a parameter set in the POWHEG BOX data file. Within this loop, gen_btilde and 
gen_uborn_idx are called in sequence. After that, radiation kinematic variables are set up 
randomly. The program then loops over all valid radiation regions (i.e. rad_kinreg values). 
For the ISR radiation region, the initial-state radiation phase space is generated, with a 
call to gen_real_phsp_isr_radO, and for the final state a call to gen_real_phsp_f sr_radO 
is performed. The task of effectively increasing the norms is performed in the routine 
inc_norms. The phase space generation routines are slight variants of the phase space 
routines previously encountered. They perform a similar task, but are dependent upon 
the rad_kinreg setting (rather than the kn_emitter value) and furthermore they compute 
the kinematics starting from the values of kn_csitilde, kn_y and kn_azi (details are 
found in the gen_real_phsp . f code). The inc_norms subroutine first sets factorization 
and renormalization scales for radiation (the set_rad_scales call), then computes the 
Born and real cross section. The real cross section is multiplied by the Jacobian J". The 
upper bounding function U TT is returned by the function pwhg_upperb_rad ( ) , and the ratio 

is formed. Its maximum gives the Nj r b (£,y) normalization. Notice that the subroutines 
sigborn_rad and sigreal_rad return respectively the value of the Born cross section for 
the current underlying Born (i.e. for the underlying Born indexed by rad_ubornidx), and 
the R TT real cross section. Thus, the subroutine sigreal_rad is similar to sigreal_btl, 
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but it only computes the cross section contributions that share the underlying Born stored 
in racLubornidx and the singular region stored in rad_kinreg. 

The first two indexes, & and yi, in the array rad_csiynorms represent the step number 
of the stepwise function 7Vj£(£,y) (i.e. they are integer functions of £ and y respectively). 
Their form can be found in the code, but may be subject to future modifications. 

For the purpose of tuning the choice of the upper bounding function, each evaluation 
of formula (7.20) is histogrammed, and the histogram is printed in TOPDRAWER format at 
the end of the upper-bound evaluation, in the file pwghistnorms .top. The efficiency in 
the generation of radiation will depend upon the shape of this histogram. It is roughly 
estimated by the ratio of the average value of formula (7.20) over its maximum. Highest 
efficiencies are achieved if the histogram goes sharply to zero near the maximum value of 
the abscissa. Lowest efficiencies are characterized by histograms with long tiny tails. 

7.1.2 The gen_radiation routine 

This routine is invoked from pwhgevent, after the call to gen_btilde and gen_uborn_idx. 
It loops through the valid radiation regions (i.e. the allowed rad_kinreg values) and it 
calls either the gen_rad_isr or the gen_rad_f sr routines, that generate and store in the 
global variables kn_csi, kn_y and kn_azi a set of kinematics radiation variables. It also 
returns, in its argument, the value of the radiation transverse momentum squared, t, that 
is defined in the function pwhg_pt2. If t is the largest generated so far, the kinematics 
radiation variables and the rad_kinreg value are saved in local variables, because they are 
the candidate for the highest bid method (discussed in appendix B of ref. [2]). At the end 
of the loop, if no call has generated any radiation, the routine exits, after setting kn_csi 
to zero, which signals the generation of a Born-like event. If radiation was generated, the 
saved values of the radiation variables are restored in global variables, and the appropriate 
phase-space generation routine is invoked. The sigreal_rad routine is invoked again, 
followed by a gen_real_idx call. Besides returning R", sig_real_rad also stores each cross 
section contribution, so that, after the radiation kinematics is generated, the corresponding 
flavour structure can also be generated with a probability proportional to each cross section 
contribution. This is what the gen_real_idx call does. The index (i.e. the air) of the 
corresponding real flavour structure is stored in the variable rad_realalr. 

The subroutine sigreal_rad, as stated earlier, is similar to sigreal_btl, but it 
only computes the cross section contributions that share the underlying Born stored in 
rad_ubornidx and the singular region stored in rad_kinreg. It also takes care to avoid 
generating gluon splittings into heavy quark pairs below threshold. More precisely, if the 
radiation k^, returned by the function pwhg_pt2(), is below rad_charmthr2 for g —> cc or 
below rad_bottomthr2 for g — > bb, the corresponding result is set to zero. 

7.1.3 The gen_rad_isr and gen_rad_f sr routines 

These routines generate a p T value according to the Sudakov form factors in eq. (7.6). 
They make essential use of the function pt2solve(pt2,i), that represents the function 
log A(p° ld )/A(p T ). The expression of logA(p T ) is as given in eq. (D.18), in the case of 
initial-state radiation, or as given in eqs. (C.6) or (C.10) (in this last case, depending upon 
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the form of the upper-bounding function used, controlled by the value of the integer variable 
iupperfsr). The value of p^ d corresponds to the last vetoed p T . The value log A(p° ld ) is 
represented by the variable xlr. Thus, solving pt2solve for zeros represents the first step 
of each iteration of the veto procedure. 

We examine now in detail the case of final-state radiation, with iupperf sr=l, that 
also illustrates how the other cases work. Looking at the function pwhg_upperb_rad, we 
see that, for iupperf sr=l, the upper bounding function (up to the normalization factor) 
has the form 

U™ c£ w ,„ , 

where Og w is the variable-flavour two loop expression for the strong coupling constant 
used for radiation generation throughout the POWHEG BOX program. It is set by a call to 
set_rad_scales (the choice of scales are discussed in detail in appendix E). In appendix C, 
it is shown how to deal with this form of the upper bounding function for the case of one 
loop, constant flavour a s , and for constant iV rr (i.e. not dependent upon £ and y). We thus 
begin by considering the upper bounding function 

ry rad 

u " = N fl ^rr^y > ( 7 - 22 ) 

with a s ad given by the one-loop expression (see eq. (C.3)) 

< d M= d / „ a • (7-23) 
^o log -ft- 

rad 

We must choose b r ad and A ra d in such a way that 

al^^a™^), (7.24) 

in all the range fi > p™ m , where p™ m is the minimum allowed p T for radiation. If this 
inequality is fulfilled, we will have U rr ^ U" in all the relevant range. The value of 6 ad is 
taken equal to 

r = ^Jr^ , (7-25) 

while A ra d is computed and stored in the global variable racLlamll by a call to the sub- 
routine init_rad_lambda at initialization stage, from the routine init_phys. 

The routine gen_rad_isr proceeds by initializing the variable unorm to the value of Ny , 
stored in the rad_norms array. The variable unorm is also made available, via a common 
block, to the function pt2solve. In the same way, the value of kt2max, appropriate to the 
current kinematics and radiation region, is computed and made available to the pt2solve 
routine, together with the value of A ra( j and the number of flavour (i.e. 5) to be used in 
&Q ad . At this stage the function pt2solve is in the condition to operate properly. Its return 
value, for iupperf sr=l, corresponds to formula (C.6). The veto loop is started, with the 
variable xlr set to the log of a uniform random number < r < 1. The zero of the 
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pt2solve function is found (using the dzero CERNLIB routine), which thus corresponds to 
a p T value that solves the equation 



log(r) - logA^ rr )(p T ) = 0. 



(7.26) 



If p% (denoted by t in the POWHEG BOX) is below the allowed minimum value, a negative 
t is returned to signal the generation of an event with no radiation (i.e. with Born-like 
kinematics). Otherwise a sequence of vetoes is applied. First, the event is accepted with a 
probability 



a 



pw (pI; 



aT d {p 2 T ) 1 (7 ' 27) 
and vetoed otherwise. After this veto is passed, the distribution of eq. (7.22) has been 
corrected for the use of a s ad , instead of the correct one, ag w . At this stage, £ is generated 
(at fixed p T ): its probability distribution is uniform in its logarithm, as can be evinced 
from eq. (C.6). The value of y is obtained by solving for y the k% = p\ definition for FSR, 
i.e. eq. (C.2). At this stage we can compute a further veto, accepting the event with a 
probability 

fbK ' (7.28) 



N fl 

which is the number returned by the subroutine uboundf ct. After this veto, £ and y have 
been generated with probability 



exp 



-J u rr e(k T -p T ) 



d£ dy deft 



2ir U rr d£dy. 



(7.29) 



At this stage, the Born cross section is computed with a call to sigborn_rad, a uniform 
azimuth for radiation is also generated and also sigreal_rad is called to compute the real 
cross section. One now vetoes again accepting the event with a probability 



J"R TT 1 



(7.30) 



and, after this, the £, y and <p variables have been generated according to the probability 

V(k T -p T ) d<J> rad 

1 h 



exp 

which is the desired result. 



-—9(k T -p T )d$ rad 



B 



(7.31) 



7.2 Remnant radiation 

Within the subroutine pwhgevent, the generation of an event with the remnant component 
of the cross section is carried out as follows. First the subroutine gen_sigremnant is 
invoked. This subroutine uses the routine gen to generate a point in the full phase space, 
distributed with a probability proportional to the sigremnant cross section, using the 
grids previously prepared, as described in section 6. The phase space point remains stored 
in the kinematics global variables. After that, the gen_remnant subroutine is invoked. 



- 39 - 



This subroutine generates the flavour structure of the current event with the appropriate 
probability. This is possible because the subroutine gen_remnant stores in the global 
arrays rad_damp_rem_arr and rad_reg_arr (defined in pwhg_rad . h) each contribution to 
the cross section for the last kinematics point computed. The genjremnant subroutine picks 
a contribution with a probability proportional to the values stored in these arrays. If the 
contribution is a remnant (described in section 5), its index is stored in rad_realalr, and 
the corresponding underlying Born and emitter is found. The radiation phase space is thus 
generated again with this value of the emitter, and the same values for the three parameters 
used to parametrize the radiation phase space in sigremnant. These parameters are stored 
by the sigremnant subroutine in the global array rad_xradremn. The recalculation of the 
radiation phase space is necessary, since only in the case when gen_remnant picks the 
last contribution computed, the phase space would already have the appropriate settings. 
The genjremnant subroutine also returns in its integer argument the value 1 for remnant 
contributions or the value 2 for regular contributions. If the contribution is from a regular 
part, its index is retrieved and stored in radjrealreg, and the ISR phase space is used, 
since this is the one we have chosen to use for all regular contributions. 

7.3 Completion of the event 

For simplicity, in the PDWHEG BOX, one always assumes that there is an azimuthal symmetry, 
so that, in the generation of the Born phase space, one can always require that some 
reference particle in the final state lies on the xz (or yz) plane, where z is the direction 
of the beam axis. At the end of the event generation, a random azimuthal rotation of the 
whole event is performed. This is done within the pwhgevent routine, through a call to 
the subroutine add_azimuth. 

Besides setting up the kinematics and the flavour structure, in order to pass the event to 
the Les Houches Interface for User Processes [38] (LHIUP from now on) , we must also decide 
up to which scale the subsequent (SMC generated) shower should start. In case of a btilde 
generated event, this scale should coincide with the radiation transverse momentum. In 
case of remnant or regular contribution, this choice is to some extent ambiguous. In order 
to maintain some continuity of the remnant events with the btilde events, we also set this 
scale to the radiation transverse momentum. For regular contributions, this value is better 
decided on the basis of the specific process, and an appropriate function pt2max_regular 
should be provided by the user, in the file pt2maxreg.f . The global variable rad_pt2max is 
set to the maximum p T for the subsequent shower. It will be used in the LHIUP interface 
to set the variable SCALUP. 

8. The Les Houches interface for user processes 

At last, the generated event is put on the LHIUP interface. The scale for subsequent 
radiation is setup, and colours are assigned to the incoming and outgoing partons. For B 
generated and remnant events, this task is carried out by the subroutine gen_leshouches. 
For regular remnants, a special routine, gen_leshouches_reg, does the job and should be 
provided by the user in the file LesHouchesreg.f . The different treatment in the two cases 
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is due to the fact that, in the case of B and remnant events, we have a standard method 
to assign colour, that is correct in the singular region. For regular contributions, instead, 
other methods should be used, like resorting, for example, to the planar limit of the cross 
section formulae. 

In general, there is much room for improvement in the technique used for colour as- 
signment [2]. We do not consider this a crucial problem at the present stage. However, if 
the need of a better colour treatment will emerge, it is clear that the user should provide 
more colour information. We thus limit ourselves, in the present case, to an approximate 
colour implementation that is general enough to be process independent. What we do is 
to assign colour on the basis of the underlying Born configuration first. Then, depending 
upon the region we are considering, we assign the colour of the real emitter and radiated 
parton as if a collinear splitting process had really taken place. In the planar limit, this 
yields a unique colour prescription for the emitter and the radiated parton, except for the 
case of a gluon splitting into two gluons, that yields two possible colour assignments with 
equal probabilities. In figures 3 and 4 two particular cases are illustrated. 



Figure 3: Colour assignment for a singular region corresponding to a quark radiating a collinear 
gluon in the final state. 



Figure 4: The two alternative colour assignments for a singular region corresponding to a gluon 
radiating a collinear gluon in the final state. 

The routine gen_leshouches begins with a call to born_lh, that sets up a few event- 
specific LHIUP parameters, like the flavour, status and mothers of the underlying Born 
incoming and outgoing particles. In order to understand this part of the code, the reader 
should refer to the specifications of the LHIUP [38] . The born_lh subroutine sets up also 
the colours of the underlying Born process, by calling the borncolour_lh subroutine. This 
subroutine is process dependent, and must be provided by the user. It should return, in the 
LHIUP, a planar colour connection, with a probability proportional to its Born contribution 
in the planar limit. For the most simple processes, like for example Z production, there 
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is one single planar colour connection, i.e. the incoming quark and antiquark should have 
complementary colours. For more complicated processes, more colour structures can arise. 
In the planar limit, different colour structures do not interfere, so it is possible to generate 
a single colour structure with a probability equal to the corresponding contribution to 
the cross section. In some cases, it may be useful to include also some colour-suppressed 
contributions. This is the case, for example, in heavy flavour production, where the leading 
colour configuration always leads to a heavy-flavour pair in an octet colour state. Singlet 
production may be more suited to the direct production of bound states, and so, it may 
be better to include it. This can be done, as long as different colour configurations do not 
interfere with each other. 

In the case when £ = 0, a Born event is produced, with a very low value of the SCALUP 
variable, so that no further radiation is generated by the SMC. Within the gen_leshouches 
subroutine, this is achieved by calling born_lh, by copying the kn_pborn momenta on the 
LHIUP (a task carried out by the subroutine moment a_lh) and by setting SCALUP to the 
minimum radiation transverse momentum. If £ / 0, radiation has taken place. The routine 
born_lh is called first, and one more parton is added with the flavour of the radiated parton; 
the leg corresponding to the emitter in the real graph is assigned the correct flavour. The 
colours are assigned on the basis of the Born colour already stored in the LHIUP. The 
subroutine set colour jrad is used to perform this task for both initial- and final-state 
radiation. 

A final task of the LHIUP subroutine is to put on the interface the intermediate 
resonances. The LHIUP specifies how resonances should be put in the interface. This 
information should be made available to the shower program, since resonance masses must 
be preserved by the shower. This is achieved by calling the user routine resonances JLh, 
which calls the add_resonance routine for each particle id of intermediate resonances, 
specifying also its decay products. 

9. Conclusions 

In this work, we have introduced and documented the POWHEG BOX, a computer framework 
for the construction of a POWHEG implementation of any given NLO process. The POWHEG 
BOX code is available at the http://moby.mib.infn.it/~nason/POWHEG/. 

In the POWHEG BOX package, the user can found three directories: W, Z and VBF_H. They 
contain the code for W, Z and Higgs boson production in vector boson fusion, and can 
serve as template for any further process that a user may want to implement. 

We would like to emphasize that the POWHEG BOX is a tool to develop programs. It is 
not something ready to run out of the box. Thus, even in order to compile the examples, 
the user should examine the Makefile, and make sure that the pdf and jet libraries are 
available in the system and are linked with the correct path. The copy of the code present 
in the repository is the SVN current version, marked with its version number. From time 
to time, we will make new SVN versions available as we implement new processes ourselves, 
or following important code improvements. 
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A byproduct of our work is the implementation of the NLO corrections for an arbi- 
trary hadronic collision process, within the Frixione, Kunszt and Signer (FKS) subtraction 
scheme. The authors of ref. [32] have also proposed a general FKS implementation for 
e + e~ processes, and they are currently extending it to hadronic collisions [39]. We stress 
that in our approach, however, the role of the NLO calculation is only meant to test the 
consistency of the implementation, and no effort was made to improve its efficiency. 
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A. Soft contributions 

In this appendix we document the calculation of the soft contribution in the FKS subtrac- 
tion framework. The real cross section in the soft approximation is given by 



2c 



1Z = 47ra s /z R < 



En ki " k i k V k i n 



+ K*, (A.l) 



where Bij is the colour correlated Born cross section of eq. (2.6), and Cj is the Casimir 
invariant for the i th leg. VJ has no singularities as the momentum of the radiated parton, 
k, goes to zero. 

A.l Soft phase space 

The phase space in the soft limit always factorizes as 

2fc (27r) d - 1 

We write now 

d^k = dki dk 2 d d ~ 3 k ± = dki dk 2 dk L k~ 2e n l ~ 2e , (A.3) 
where we have set d = 4 — 2e, and £l a is the solid angle in a dimension 

-r(l + a/2)~ ^T{a) => u ~ l r(l - 2e) ' [ ' 

Turning eq. (A.3) into polar coordinates we get 

' "' ; fc^ 2e (sin^sinc/»)- 2e ciA;ocicos0dc/», (A.5) 



2fc°(2 7 r) d - 1 r(l-2e) (2vr) 3 
where we have defined 

ki = kocos6, k 2 = ko sin 6 cos (f>, k± = ko sin 9 sin <\>. (A.6) 
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Since > 0, this means that < <fi < ir, and that only even quantities can be integrated 
in this way. In other words, k± should not be confused with k^ (k% is no longer available 
at this stage). Inserting 



2 ' 



(A.7) 



this becomes 
d^k 



2A; (2vr) 



d-l 



(4^(1-6) 
T(l - 2e) 



1 s 



£ 1_2e (sin 8 sin <^)" 2e ^ d cos 6 d<f> . (A.8) 



>-2e 



(2vr)3 4 

By writing the 7£ term as £~ 2 (£ 2 7£), we notice that (£ 2 7£) has a finite limit as £ — ?■ 0. The 
£ integration is performed by separating first 



-l-2e 



-2c 



2e 
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log_£ N 
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where the <5(£) term yields the soft contribution. We thus have that the integral of the 
soft-divergent part of 1Z is given by 
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(4tt)T (1 - e) 
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(2vr) s 



d cos dd>( sin 8 sin <£) 2e 
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(A.10) 



where 7£ s is now independent upon £ (the dependence on £ of A; is canceled by the £ 2 term 
in the numerator). Collecting the normalization factor of eq. (2.12) in front, we get 

2e J 2vr 
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We now define 
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so that 
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(A.13) 
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We introduce then our basic integral 



I{jP-iO) = I dcosd — (sin#sin< 

J 7T 



-2e 



s£ 2 p 



^d(p, ?) + hip, q) + chip, q)- 



4 p ■ kq ■ k 

(A.15) 

The expression of I(p, q) will be substantially different for the case when both p and q are 
massless, when one is massive and one massless, and when both are massive. 
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A. 2 One massless and one massive particle 

Consider two momenta p and to, with p 2 = and m 2 7^ 0. We want to evaluate I(p,m). 
We separate out the collinear divergent component from the eikonal factor 
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p ■ n 
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(A.16) 
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where the square bracket term has no collinear singularities. Assuming n along the time 
direction, we have: 



n ■ p 



1 



and 



4 p ■ k n ■ k 1 — cos 9 ' 
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1 — cos 8 e 
so that 

Id(p,m) = -1. 

The remaining integral has no collinear singularities so that we can write 
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we have 
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(A.17) 
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Thus, eq. (A. 12), in the case where k\ is massless and k2 is not, is given by 



•-12 



7T 



Q_ 

si 2 



2\ e 



— I / d cos — (sin sin <p) 2e 
2eJ J vr v v> 



s£ 2 ki ■ k 2 
4 k\ ■ k k 2 ■ k 



Q 2 



A B r, 

-2+- + C, 



Q 2 



TT- 



1 + elog — r + - log —r - — e 



si 



Te )I(h,k 2 ) 



where 



A = 
B = 

C = 



1 Q 2 



1 



Io(h,k 2 ) , 



2^ 

si 2 



7T 

IF 



1 Q 2 1 

- 2 J o( fc i, fo) log ^2 ~ 2 7e ^ 1 ' k2) ■ 



(A.25) 

(A.26) 
(A.27) 

(A.28) 



-45 - 



A. 3 Two massless particles 

Using the identity 

h-k 2 



h ■ (h + k 2 ) 
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k 2 ■ (h + k 2 ) 
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that holds if k\ = and k\ = 0, we can immediately obtain the expression of I(k\, k 2 ) for 
two massless momenta 



I(ki,k 2 ) = I(ki, ki + k 2 ) + I(k 2 , h + k 2 ), 
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and use the results of the previous subsection for the two terms on the right hand side. We 
can write 
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A. 4 Two massive particles 

In case both fci and fc 2 are massive, we define 
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and get (neglecting now e 2 terms) 
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with 



B = ~I (k 1 ,h 2 ), (A.39) 
C = -^Io(k 1 ,k 2 )log^-^h(k 1 ,k 2 ), (A.40) 



and 



/o(tl , fa)= l log i±|, /S^l-^M.. (A.4!) 

The expression for is defined by the equations below 
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(«i)o (fojO 

The calculation of this integral is long and cumbersome. However, its correctness can be 
easily checked numerically, once the analytic answer has been obtained. Codes for checking 
the soft integrals are included in the Notes subdirectory of the PDWHEG BOX. 
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A. 4.1 Two massive particles with equal momenta 

In the particular case when k± = ki = p, p 2 / 0, we have 

Io(p,p) = 2, 

so that 
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B. Collinear limits 

We list here the relationship between real- and Born-level squared amplitude, summed 
over the final-state and averaged over initial-state colours and spins, and divided by the 
appropriate flux factor. 

B.l Initial-state radiation 

We call p the momentum of the incoming parton and k the momentum of the parton that 
enters the underlying Born process. Thus, (p—k) is the momentum of the on-shell radiated 
parton, and k 2 < 0. We define 
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B.2 Final-state singularities 

We call k the momentum of the splitting parton, and write 

I -, 2 



We have 
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In eq. (B.9), it is assumed that p is the momentum of the gluon. 
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C. Upper bounding functions for FSR 

We use an upper bounding function of the following form (for ease of notation we write a s 
instead of a| ad ) 



U(£,y) = N 



a s (k 2 T ) 

f(i-y)' 



with 
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s being the partonic CM energy squared, and 
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The ranges of £, y and <f> are given by 
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where M rec is the mass of the system recoiling against the emitter and emitted partons 
(see eq. (5.49) or ref. [2]). We want to generate p T uniformly in 
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Trading y for fej (see eq. (C.2)), we get 



-logA^)(p T ) = 2vriV /^^ [ 
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The equation r = A^(p T ) is translated into logr = log A^ u \p T ), and solved numerically 
for p T . Once p T is generated, we generate £ uniformly in log£ (see second member of 
eq. (C.6)) within the limits 

y < ^ U , (C.7) 

and then use 

Px = ^ 2 (i-y) (C8) 

to obtain y, while is generated uniformly between and 2tt. 

More options are offered in the POWHEG BOX. They are both related to the use of the 
upper bounding function 

U&y) = N ^ ^, (C.9) 
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where p = \Jp%j s. The second option implemented in the POWHEG BOX makes use of the 
upper bounding function in eq. (C.9) multiplied by a factor of £. One then uses the same 
expression (C.10) for the generation of p T , and uses a further veto, accepting the event if a 
given random number is less than £, to implement the extra factor of £. The third option 
is similar to the second, but with an extra factor of 1 — £(1 — y)/2. 

D. Upper bounding functions for ISR 

The upper bounding function is (where x = 1 — £, so that the singular limit is reached 
when x — >■ 1) 

(l-x)(l-y 2 

with 
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Sb being the underlying Born CM energy squared. The range of U(x, y) must cover the 
range of the radiation variables for the given underlying Born configuration. A practical 
restriction for the range of U(x, y) is 

p^x^l, k 2 T ^k 2 T ^, (D.3) 

where 
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in the given range. We assume 0^0^ 2ir. Trading y for k% we find 
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where 
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The x integration can be performed to yield 
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where 



jdxjdyj^ dct>U(x,y)9(k T -p T ) = J 2 F(^) , (D.9) 



F(^) = 27 rATa s (^) log ^ + P + ^El g . (D.10) 

y/x + - p - y/xZ - p 



We observe that 
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The dk 2 integral of V can be performed analytically, yielding 
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One generates p T uniformly in 

A(p T ) = exp 



fc Tmax (7i- 2 _ 
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and then use the veto method to get the p T distributed according to A^^t). 

The following variant of this procedure has been introduced in ref. [4]. We have used 
the bound 
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To improve the behaviour for small x effects it may be convenient to use instead 

U(x,y) = N a »ffi 
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as upper bounding function. As derived in ref. [3] 
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In this case too V(k%) satisfies a simple upper bound 

V(kl) < irNa s (kl) log ^ , (D.23) 
that can be used for fast generation by vetoing. 

E. Choice of scales 

E.l Scales and couplings for the inclusive cross section 

In the evaluation of B and of the remnant function, a user provided, process dependent sub- 
routine set_f ac_ren_scales (muf ,mur) sets the factorization and renormalization scales. 
It should only depend upon the underlying Born kinematics (thus, for example, it can- 
not depend upon the radiation transverse momentum). It is called by the POWHEG BOX 
subroutine setscalesbtilde (called during the evaluation of the B function and of the 
remnant cross section) that stores the square of the factorization scale, the square of the 
renormalization scale and the strong coupling constant in three global variables of the st_ 
common block: st_muren2, st_mufact2 and st_alpha. By inspecting the subroutine we 
see that the factorization and renormalization scales are set equal to the square of the 
scales returned by the set_f ac_ren_scal.es subroutine multiplied by renormalization and 
factorization scale factors st_facfact and st_renfact. These factors are in turn read 
from the variables f acscf act and renscf act in the POWHEG data file, normally during the 
execution of the user initialization routine init_phys. If they are not present in the data 
file, they are set by default to 1. The function pwhg_alphas(mu2,Lambda5,n) returns the 
strong coupling constant with n light flavours, as a function of the square of the scale and 

(5) 

of A^. If called with n < it uses, as number of flavours, the number of quarks with 
mass less than the input scale. In the evaluation of the B function, the number of flavours 
st_nlight is used. This is set by the user in the init_couplings routine. 

E.2 Scales and couplings for radiation 

The choice of scales and couplings in the generation of radiation requires particular atten- 
tion, since the shape of the Sudakov peak, in the radiation transverse momentum, is deeply 
affected by it. It is discussed in detail in ref. [2]. Here we report how this is actually done 
in the POWHEG BOX. The relevant code is in the subroutine set_rad_scales (ptsq) . It is 
typically invoked with the radiation transverse momentum as argument. With this choice 
one can achieve, in some cases, complete next-to-leading logarithmic (NLL) accuracy in 
the POWHEG Sudakov form factor [3, 2]. By inspecting the code, we can see that it sets 
st_mufact2 and st_muren2 to ptsq. The factorization and renormalization scale factors, 
st_facfact and st_renfact, are not used in this context. The program also takes care 
that the factorization scale never goes below the minimum allowed value in the pdf 's. The 
strong coupling is then evaluated, the number of flavours is taken as the number of quarks 
with mass below the ptsq scale. Furthermore, the strong coupling constant is multiplied 
by the factor given in formula (4.32) in ref. [2]. This factor improves the NLL accuracy of 
the Sudakov form factor. 
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During the generation of radiation, we need a simplified one loop expression for the 
running coupling that is an upper bound of the running coupling that we use. We choose 



d _ 1 (5) _ 33-2x5 

s ^' ~ 12^ ' ( } 

'Mi 

and we fix An by requiring 

rad/,,0\ „PW/„0\ /T7> o\ 

"s (M ) = a s ) , ( E - 2 ) 
where the scale /io is the minimum allowed value for the renormalization scale, and is taken 

(5) 

equal to 2A^. The value of An is stored in the global variable racLlamll. 



F. Miscellaneous features of the code 

F.l Checking the soft, collinear and soft-collinear limits 

In the POWHEG BOX, the routine checklims, through a call to the subroutines checksoft 
and checkcoll, allows the user to check the limiting behaviours of the real squared am- 
plitudes, against their soft and collinear approximations. This routine has been used as a 
debug feature in the developing stages of the implementation of specific processes. If acti- 
vated by the flags dbg_sof ttest and dbg_colltest in the init_phys.f file, it can be used 
now to check if the Born, the colour-correlated and spin-correlated Born amplitudes (used 
to build the limiting expressions of the real amplitudes) are consistent with the real con- 
tributions, computed in the kinematic configurations where a gluon becomes soft or when 
it becomes collinear to another parton or when two quarks of opposite flavours become 
collinear. 

The double soft-collinear and collinear-soft limits are also tested. They do not depend 
upon the real squared amplitudes, but only upon the Born amplitudes. They have been 
used initially to check the consistency of the POWHEG BOX code, but they can also be used 
to perform some checks of the colour correlated Born amplitudes during the development 
of new processes. 

F.2 Names 

In the present work, we have made little references to the names of the fortran files where 
the various subroutines are stored. We assume that the reader can find them out by using 
standard command-line tools. However, while we have not spent much energy in deciding 
how to organize fortran files, the include files are rather well organized. For example, the 
include file pwhg_f 1st .h declares common-block variables that refer to flavour structures. 
All these variables start with a prefix f lst_. Other important global variables are the 
kinematics ones (kn_), those involving the scales and the strong coupling (st_), those 
involving the generation of radiation (rad_), and so on. In this way, when finding such 
variables in the code, by inspecting the include files, the reader can better follow what is 
their purpose. 
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F.3 Input variables 

The POWHEG BOX gets its input data from the routine powheginput. Upon its first in- 
vocation, this routine looks for a file named powheg. input. If it does not find it, asks 
interactively to enter a prefix, and then looks for the file "prefix "-powheg. input. It 
reads all the input file at once. Then, if invoked in the form powheginput ("string") it 
returns the (real) value associated to "string" in the input file. If no matching string is 
found, it prints a message and aborts the program. If invoked in the form "#string", in 
case no matching string is found, it returns a very unlikely value — 10 6 . This last mechanism 
is used in the POWHEG BOX to set default values. 

F.4 User files 

The files that contain the user routines are organized in the same way in all the examples 
provided with the code of the POWHEG BOX. They are: 

• nlegborn.h contains the number of legs of the Born process, nlegborn; 

• init_processes . f contains the subroutine init_processes, whose major task is to 
fill the list f lst_born and f 1st jreal; 

• init_couplings . f contains the subroutine init_couplings, that initializes process- 
dependent couplings; 

• in PhysPars.h, there is a collection of physical variables that are in common with 
many subroutines (masses, electroweak couplings, widths.. .); 

• Born . f contains the setborn subroutine; 

• in Born_phsp . f the user can found the born_phsp subroutine; 

• in real.f, the setreal routine; 

• and finally, in virtual. f, the setvirtual routine. 

A template for the analysis subroutine can be found in pwhg_analysis . f . 
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